C_____________________________________________________________________ C | C ***** **** ***** * * ****** * * | C * * * / * * * * * | C * * * \****\ ****** **** * * | C * * * * * * * * * | C ***** **** *****/ * * ****** ****** ****** | C | C Program Name: COSHELL | C | C Purpose: Linear stress, linearized buckling, natural | C vibration, nonlinear and post-buckling analysis | C of stiffened composite cylindrical shells using | C the Layerwise Shell Theory. | C | C Coded by: Samuel Kinde Kassegne | C | C Virginia Polytechnic Institute and State Univ. | C Blacksburg, VA 24061 | C September 1992. | C_____________________________________________________________________| C C_____________________________________________________________________ C C List of variables, arrays and matrices follow: 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 NMX -- Size of global stiffness matrix, used only for postbuckling C analysis, (i.e. NMX=NUX) otherwise, NMX = 10 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 DUT(NMX) -- Tangential displacements (postbuckling analysis) C DUR(NMX) -- Residual displacements (postbuckling analysis) C DU1(NMX) -- Updated(1) displacements (postbuckling analysis C DU2(NMX) -- Updated(2) displacements (postbuckling analysis) C DELU(NMX) -- Increment. displacements (postbuckling analysis) C FBAR(NMX) -- Global unblanced force (postbuckling analysis) C GUPLS(NMX) -- Converged displacements from previous iteration C for a given load step -used to trace back the C path with different arc length if negative roots C encountered (postbuckling analysis) C NOD,NODR,NODS -- Connectivity array for shell,ring and stringer C elements respectively C THETA,THETR -- Array for the ply orientations of shell,ring and C THETS stringer laminates C THICK,THICKR -- Array for the ply thicknesses of shell,ring and C THICKS stringer laminates C Z, ZZ -- Array for depth-coordinates of laminates of C shell and stiffeners C RADIUS(15) -- Array for the radius of each ply in the shell C VBSF,IBSF -- Arrays for natural (force) boundary conditions. C VBDY,IBDY -- Arrays for esst. (geometry) boundary conditions. C ILOCS,ILOCR -- Array for the location of each stringer and ring C = 0 for internal, C = 1 for external stiffener. C C ___________________ C C ii) control variables C C NEIGEN -> = 1 for eigenvalue analysis C IPROB -> = 0 for natural vibration, C = 1 for buckling C IMPFCT -> = 1 for imperfection sensitivity analysis C ISINE -> = 0 for uniform load. C = 1 for double sinusoidal load, C = 2 for single sinusoidal load (x-dir) C = 3 for single sinusoidal load (y-dir) C ISMEAR -> = 0 for unstiffened case, C = 1 for full discrete model, C = 2 for equivalent discrete model and C = 3 for smeared approach C ISTART -> = 0 for no restart, C = 1 for write only, C = 2 for read only and C = 3 for read/write options C C IPOST -> = 1 for postbuckling analysis 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 shell skin elements C = 2 for printing stresses and force resultants C of both the shell and stiffener elements C C C ___________________ C C C iii) data variables C C NELEM -- Number of shell 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 analy 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 DLAM -- Maximum arc length (postbuckling analysis) C C_______________________________________________________________________ IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER (NMX=10,NEX=10,NUX=486,NGX=150,NBX=126) COMPLEX*16 ALPHA(NEX) COMMON /SEL/ STIF(NGX,NGX),ELP(NGX),ELXY(9,2),ELFC(NGX) COMMON /GLOBAL/GSTIF(NUX,NBX),GF(NUX),W(NUX),EPL(NUX) COMMON /RIK/ DUT(NMX),DUR(NMX),DU1(NMX),DU2(NMX),GUPLS(NMX) COMMON /ANGLE/ THETS(20),WIDTH(20),THICKS(20),THETR(20),THICKR(20) COMMON /MSR/ NODR(50,4),NODS(50,4) COMMON /BCS/ VBDY(525),IBDY(525) COMMON /STIFF/ E1S,E2S,G12S,ANU12S COMMON /STRING/AS,ETS,EJS,ES,GS,BS,E3S COMMON /RING/ AR,ETR,EJR,ER,GR,BR,E3R COMMON /ELAST/ THETA(20),RADIUS(15),ELX(3) COMMON /DLDA/ DLAMDA(15),DLAM COMMON /TIT/ TITLE(20),VBSF(40) COMMON /IBS/ IBSF(40) COMMON /DEPTH/ Z(15),ZZ(15),THICK(15),SM(NGX,NGX) COMMON /EDWR/ EDWDX(25,20),EDWDY(25,20) COMMON /MSH/ NOD(50,9),X(125),Y(125),DX(15),DY(15) COMMON /CCC/ E1,E2,E3,G12,G13,G23,ANU12,ANU13,ANU23 COMMON /SMAS/ GMASS(NEX,NEX),EGNVAL(NEX) COMMON /PBUK/ DELU(NMX),FBAR(NMX) COMMON /SMEAR/ EIS,SS,ZS,EIR,SR,ZR COMMON /LOCATE/ILOCS(30), ILOCR(30),ISCOD COMMON /SPAND/ SPANX,SPANY COMMON /PAR/ NRED,ILOD,IPROB,IPOST,NIT,LL,IAXIAL,ILATER,ISHEAR COMMON /PAR1/ NEIGEN,IFLAG,NPLY,NPLS,NPLR COMMON /PAR2/ IMPFCT,ISINE,ISMEAR DIMENSION T(6) C COMMON /WORKSP/ RWKSP C REAL RWKSP(10) C CALL IWKIN(10) CHARACTER*14 FLINP,FLOUT DATA IN,IT /5,6/ DATA NRMAX,NCMAX/NUX,NBX/ C______________________________________________________________________ C C.....READ IN THE FINITE ELEMENT MESH DATA C______________________________________________________________________ WRITE(6,222) WRITE(6,325) WRITE(6,337) 337 FORMAT(//,' Enter COSHELL Input File --->') READ (*,334) FLINP 334 FORMAT(\A) WRITE(*,335) 335 FORMAT(/,' Enter COSHELL Output File -->') READ (*,336) FLOUT 336 FORMAT(\A) OPEN(UNIT=5,FILE = FLINP,STATUS='OLD') OPEN(UNIT=6,FILE = FLOUT,STATUS='UNKNOWN') OPEN(UNIT=7,FILE = 'matrix.out',STATUS='UNKNOWN') OPEN(UNIT=8,FILE = 'DISP.out',STATUS='UNKNOWN') WRITE(6,222) READ(IN,300) TITLE READ(IN,*) NELEM,NPE,NNODE,MAXIT,NLOAD,NEIGEN,IPROB,IPOST,IR READ(IN,*) PX,PW,EPSLON,NPRT,NRED,ISTRES,ICODE,ILOD,DLAM,ISTART READ(IN,*) IAXIAL,ILATER,ISHEAR,IREAD,KNOD,ISMEAR,IMPFCT,ISINE DO 20 N = 1,NELEM 20 READ(IN,*) (NOD(N,I), I=1,NPE) READ(IN,*) (X(I),Y(I),I=1,NNODE) READ(IN,*) NBDY,NBSF IF(NBDY.NE.0.AND.IREAD.EQ.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 C.... READ IN SHELL SKIN LAMINATE DATA C_____________________________________________________________________ IF(ISMEAR.EQ.1) THEN READ (IN,*) NPLY,NPLS,NPLR,ISCOD ELSE READ (IN,*) NPLY END IF READ (IN,*) (THETA(I), I=1,NPLY) READ (IN,*) (THICK(I), I=1,NPLY) TH = 0.0 DO 24 I = 1,NPLY 24 TH = TH + THICK(I) NIF = NPLY+1 IF(ISCOD.LE.1) NPLR = 0 IF(ISMEAR.EQ.1) NIF = NPLY+1+NPLS+NPLR IF(IMPFCT.EQ.0) THEN NDF = 3*NIF ELSE NDF = 4*NIF END IF JC =1 IF(IREAD.EQ.1) THEN DO 45 K = 1,KNOD READ(5,*) NODE,IC J = (NODE-1)*NDF IF(IC.EQ.111.OR.IC.EQ.101.OR.IC.EQ.110.OR.IC.EQ.100) THEN DO 15 I = 1,NIF IBDY(JC) = J+I VBDY(JC) = 0.0 JC = JC +1 15 CONTINUE END IF IF(IC.EQ.111.OR.IC.EQ.110.OR.IC.EQ.011.OR.IC.EQ.010) THEN DO 25 I = 1,NIF IBDY(JC) = J+NIF+I VBDY(JC) = 0.0 JC = JC +1 25 CONTINUE END IF IF(IC.EQ.111.OR.IC.EQ.101.OR.IC.EQ.011.OR.IC.EQ.001) THEN DO 35 I = 1,NIF IBDY(JC) = J+2*NIF+I VBDY(JC) = 0.0 JC = JC +1 35 CONTINUE END IF IF(IC.EQ.222) THEN DO 37 I = 1,NIF IBDY(JC) = J+I VBDY(JC) = 1.0 JC = JC +1 37 CONTINUE END IF 45 CONTINUE END IF C_______________________________________________________________________ C C ENTER GEOMETRICAL IMPERFECTIONS OF THE SHELL AS DATA. C (IF IMPFCT = 1) C_______________________________________________________________________ IF(IMPFCT.EQ.1) THEN READ (IN,*) GEOMAG IF(IREAD.EQ.0) JC = NBDY+1 DO 85 I = 1,NNODE J = (I-1)*NDF DO 65 K = 1,NIF IBDY(JC) = J + 3*NIF + K VBDY(JC) = GEOMAG JC = JC + 1 65 CONTINUE 85 CONTINUE NBDY = JC-1 END IF READ (IN,*) E1, E2, E3, G12, G13, G23, ANU12, ANU13,ANU23,RHO C_______________________________________________________________________ C C.... READ IN STIFFENER DATA C C ISMEAR = 1 ; USE FULL DISCRETE APPROACH. C ISMEAR = 2 ; USE DISCRETE APPROACH C ISMEAR = 3 ; USE SMEARING APPROACH C C_______________________________________________________________________ IF(NPE.EQ.4) THEN NPB = 2 ELSE NPB = 3 END IF READ (IN,*) NSTRG, NRING, SPANX,SPANY,(RADIUS(I),I=1,NIF) IF(ISMEAR.EQ.0) GO TO 535 IF (ISMEAR.EQ.1) THEN READ (IN,*) (ILOCS(I), I=1,NSTRG), (ILOCR(I),I=1,NRING) READ (IN,*) (THETS(J), J=1,NPLS), (THETR(J),J=1,NPLR) READ (IN,*) (THICKS(J),J=1,NPLS), (THICKR(J),J=1,NPLR) DO 38 N = 1,NSTRG 38 READ (IN,*) (NODS(N,I),I = 1,NPB) DO 40 N = 1,NRING 40 READ (IN,*) (NODR(N,I),I = 1,NPB) END IF IF(ISMEAR.EQ.2) THEN READ (IN,*) AS,DS,NPLYS,(ILOCS(I),I= 1,NSTRG), * AR,DR,NPLYR,(ILOCR(I),I= 1,NRING) C.....STRINGER DATA FOLLOWS READ (IN,*) (THETS(J),J=1,NPLYS) READ (IN,*) (WIDTH(J),J=1,NPLYS) READ (IN,*) (THICKS(J),J=1,NPLYS) READ (IN,*) E1S,E2S,G12S,ANU12S BS = WIDTH(1) DO 32 N = 1,NSTRG 32 READ (IN,*) (NODS(N,I),I = 1,NPB+1) CALL COLAM (NPLYS,ES,GS,ETS,E3S,THETS,WIDTH,THICKS) C.... RING DATA FOLLOWS READ (IN,*) (THETR(J), J=1,NPLYR) READ (IN,*) (WIDTH(J), J=1,NPLYR) READ (IN,*) (THICKR(J),J=1,NPLYR) READ (IN,*) E1S,E2S,G12S,ANU12S BR = WIDTH(1) DO 42 N = 1,NRING 42 READ (IN,*) (NODR(N,I),I=1,NPB+1) CALL COLAM(NPLYR,ER,GR,ETR,E3R,THETR,WIDTH,THICKR) END IF IF (ISMEAR.EQ.3) THEN READ (IN,*) ES,GS,EIS,AS,EJS,SS,ZS,RHOS,ILCS READ (IN,*) ER,GR,EIR,AR,EJR,SR,ZR,RHOR,ILCR END IF C______________________________________________________________________ C C.... COMPUTE THE HALF BAND-WIDTH C_____________________________________________________________________ 535 NEQ = NNODE*NDF NN = NPE*NDF 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 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,SPANX,RADIUS(NIF) WRITE(IT,362) E1,E2,E3,G12, G13, G23,ANU12,ANU13,ANU23 IF(NEIGEN.EQ.1) THEN IF(IPROB.EQ.1) THEN LCODE = 5 ELSE LCODE = 4 END IF END IF IF(IPOST.EQ.1) LCODE = 3 IF(IMPFCT.EQ.1) LCODE = 6 IF(IPOST.NE.1.AND.NEIGEN.NE.1.AND.IMPFCT.NE.1) THEN IF(EPSLON.EQ.1.0) THEN LCODE = 1 ELSE LCODE = 2 END IF END IF IF(NBSF.NE.0) ISINE = 4 WRITE(IT,353) LCODE IF(LCODE.EQ.3) WRITE(IT,343) IR IF(LCODE.EQ.3) WRITE(IT,333) ISTART WRITE(IT,363) ISINE WRITE(IT,273) ILOD WRITE(IT,348) NRED,ICODE,ISTRES,MAXIT,EPSLON IF(NPRT.NE.1) GO TO 344 WRITE(IT,315) WRITE(IT,605) (I,X(I),Y(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 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 C.... WRITE OUT STIFFENER DATA C_______________________________________________________________________ IF(ISMEAR.EQ.1) THEN WRITE(IT,463) WRITE(IT,592) NSTRG,NRING,NPLS,NPLR,(THICKS(I),I=1,1), + (THICKR(I),I=1,1),E1,E3,G13 WRITE(IT,593) (ILOCS(I),I=1,NSTRG) WRITE(IT,594) (ILOCR(I),I=1,NRING) ELSE IF(ISMEAR.EQ.2) THEN WRITE(IT,461) WRITE(IT,392) NSTRG,NRING,ILOCS(1),ILOCR(1), + ES,GS,E3S,AS,ER,GR,E3R,AR ELSE IF(ISMEAR.EQ.3) THEN WRITE(IT,462) WRITE(IT,492) SS,ZS,RHOS,ILCS,ES,GS,AS,EIS,EJS, + SR,ZR,RHOR,ILCR,ER,GR,AR,EIR,EJR END IF C_______________________________________________________________________ C C.... PROCESS INPUT DATA C_______________________________________________________________________ C..... CALCULATE THE THROUGH-THICKNESS STIFFNESSES,[A],[D],[DB],[F] ETC C..... (VALID FOR ALL SECTIONS AND ITERATIONS) NPE1 = 2 C______________________________________________________________________ CALL SECOND(T(1)) 635 CALL STFPLY(NPLY,DEP,Z,THICK,THETA,0) CALL SECOND(T(2)) C______________________________________________________________________ IF(IPOST.NE.1.OR.ISTART.LT.2) THEN DO 95 IJ = 1,NUX IF(IPOST.EQ.1) THEN GUPLS(IJ) = 0.0 DELU(IJ) = 0.0 END IF 95 W(IJ) = 0.0 DLAMDA(1) = 0.0 DARC = 0.0 ITERD = 1 ITERP = 1 ICNVRG = 0 IERR = 0 END IF C_______________________________________________________________________ C C.... READ PARAMETERS FOR A RE-START ANALYSIS C (Only for nonlinear nalysis and ISTART.GE.2) C C_______________________________________________________________________ IF(IPOST.EQ.1.AND.ISTART.GE.2) THEN REWIND 10 READ (10,101) ITERP,DARC,DLAM,DLAMDA(1) READ (10,102) (W(I),I = 1,NEQ) DO 117 I = 1,NEQ 117 GUPLS(I) = W(I) 101 FORMAT(I3,3(E10.3,1X)) 102 FORMAT(6(E10.3,1X)) END IF C_______________________________________________________________________ C C.... ITERATIONS FOR NON LINEAR OR POST-BUCKLING ANALYSIS C_______________________________________________________________________ DO 55 LL =1,NLOAD NITER = 0 IF(NEIGEN.EQ.1) GO TO 158 IF(IPOST.EQ.1) THEN DO 82 I = 1,NEQ 82 DELU(I) = 0.0 END IF 158 NITER = NITER+1 NIT = NITER IF(NITER.LT.MAXIT) GO TO 159 CALL SOVER GO TO 94 159 CONTINUE C_______________________________________________________________________ C C.... INITIALIZE GLOBAL MATRICES C_______________________________________________________________________ DO 72 I= 1,NUX GF(I) = 0.0 IF(IPOST.EQ.1) FBAR(I)= 0.0 DO 72 J=1,NBX 72 GSTIF(I,J) = 0.0 IF(NEIGEN.EQ.1) THEN DO 74 I = 1,NEX DO 74 J = 1,NEX 74 GMASS(I,J) = 0.0 END IF IF(ISMEAR.EQ.1.AND.NITER.GT.1) THEN CALL STFPLY(NPLY,DEP,Z,THICK,THETA,0) END IF C_______________________________________________________________________ C C.... BEGIN DO LOOP ON NUMBER OF ELEMENTS C_______________________________________________________________________ DO 150 N=1,NELEM NDF = 3*(NPLY+1) NN = NPE*NDF L = 0 DO 100 I = 1,NPE NI=NOD(N,I) ELXY(I,1) = X(NI) ELXY(I,2) = Y(NI) LI = (NI-1)*NDF DO 100 J = 1,NDF LI = LI+1 L = L+1 ELFC(L) = W(LI) 100 CONTINUE IF(ISMEAR.EQ.1) THEN NIS = (NPLY+NPLR+NPLS+1) NI = NPLR IF(ISCOD.EQ.0) NI = 0 DO 28 I = 1,NPE NJ = NOD(N,I) DO 28 J = 1,3 DO 28 K = 1,NPLY+1 II = (NJ-1)*NIS*3 + (J-1)*NIS + K + NI IK = (I-1)*(NPLY+1)*3 + (J-1)*(NPLY+1) + K 28 ELFC(IK) = W(II) END IF C_______________________________________________________________________ C C.... EVALUATE THE TANGENT STIFFNESS MATRIX OF THE SHELL C C_______________________________________________________________________ C_______________________________________________________________________ CALL STFSTR(PX,PW,NPE,N,RHO,DEP,NDF,NN) C_______________________________________________________________________ IF(N.EQ.1.AND.NPRT.EQ.2.AND.NITER.LE.2) THEN DO 777 I=1,NN 777 WRITE(IT,878)(STIF(I,JJ),JJ=1,8) WRITE(IT,177) WRITE(IT,177) END IF 177 FORMAT(//,4X,72('_'),///) 878 FORMAT(8(E9.2,1X)) 873 FORMAT(6(E9.2,1X)) C_______________________________________________________________________ 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 IF(IPOST.EQ.1) THEN FBAR(NR) = FBAR(NR) + ELP(L) IF(NITER.EQ.1.OR.IR.EQ.0) THEN GF(NR) = GF(NR) + EPL(L) END IF END IF IF(IPOST.EQ.0) THEN GF(NR) = GF(NR) + ELP(L) END IF 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 C_______________________________________________________________________ C C..... CALCULATE THE ADDITIONAL STIFFNESS OF THE STRINGERS AND RINGS C C_______________________________________________________________________ IF(ISMEAR.EQ.3) THEN C_______________________________________________________________________ C C.... SUBROUTINE STFSMR USES THE SMEARING APPROACH TO MODEL STIFFENERS C_______________________________________________________________________ C IFLAG = 1 CALL STFSMR(NPLY,NPE,N,RHOS,DEP,0,IFLAG,NEIGEN,NDF) CALL STFSMR(NPLY,NPE,N,RHOS,DEP,1,IFLAG,NEIGEN,NDF) IFLAG = 2 IF(EIR.NE.0.AND.AR.NE.0) THEN CALL STFSMR(NPLY,NPE,N,RHOR,DEP,ILCR,IFLAG,NEIGEN,NDF) END IF END IF 150 CONTINUE CALL SECOND(T(3)) C_______________________________________________________________________ C C.... SUBROUTINE STFNER USES THE DISCRETE APPROACH TO MODEL STIFFENERS C_______________________________________________________________________ C IF(ISMEAR.EQ.2) THEN IFLAG = 1 ELEN = DS CALL STFNER(IFLAG,NIF,NPB,ELEN,NSTRG,RHO,NEIGEN,NDF,ILOCS) IFLAG = 2 ELEN = DR CALL STFNER(IFLAG,NIF,NPB,ELEN,NRING,RHO,NEIGEN,NDF,ILOCR) END IF C_______________________________________________________________________ C C.... SUBROUTINE STFFUL USES FULL DISCRETE APPROACH TO MODEL STIFFENERS C_______________________________________________________________________ C IF(ISMEAR.EQ.1) THEN CALL STFPLY(NPLS,DES,ZZ,THICKS,THETS,1) NIS = (NPLY+NPLR+NPLS+1) IFLAG = 1 RAD = 1E12 DO 715 N = 1,NSTRG L = 0 DO 720 I = 1,NPB NI = NODS(N,I) ELX(I) = X(NI) LI = (NI-1)*NDF IF(ILOCS(N).EQ.0) THEN DO 723 J = 1,(NPLS+1)*2 LI = LI+1 L = L+1 IF(J.LE.NPLS+1) THEN ELFC(L) = W(LI) ELSE ELFC(L) = W(LI+NPLS+NIS+1) END IF 723 CONTINUE ELSE DO 724 J = 1,(NPLS+1)*2 LI = LI+1 L = L+1 IF(J.LE.NPLS+1) THEN ELFC(L) = W(LI+NPLR+NPLY) ELSE ELFC(L) = W(LI+NIS+NPLR+NPLY+NPLR) END IF 724 CONTINUE END IF 720 CONTINUE CALL STFFUL(NPB,RHO,NDF,RAD,NRED,IPROB,N,ILOCS,DES) 715 CONTINUE CALL STFPLY(NPLR,DER,ZZ,THICKR,THETR,1) IFLAG = 2 DO 735 N = 1,NRING IF (ILOCR(N).EQ.0) THEN RAD = RADIUS(1) ELSE RAD = RADIUS(NPLY+1) END IF L = 0 DO 740 I = 1,NPB NI = NODR(N,I) ELX(I) = Y(NI) LI = (NI-1)*NDF+NIS IF(ILOCR(N).EQ.0) THEN DO 743 J = 1,(NPLR+1)*2 LI = LI+1 L = L+1 IF(J.LE.NPLR+1) THEN ELFC(L) = W(LI) ELSE ELFC(L) = W(LI+NPLS+1) END IF 743 CONTINUE ELSE DO 744 J = 1,(NPLR+1)*2 LI = LI+1 L = L+1 ELFC(L) = W(LI+NPLY+NPLR) 744 CONTINUE END IF 740 CONTINUE CALL STFFUL(NPB,RHO,NDF,RAD,NRED,IPROB,N,ILOCR,DER) 735 CONTINUE END IF CALL SECOND(T(4)) C_______________________________________________________________________ C C CARRY OUT EIGENVALUE ANALYSIS FOR NEIGEN.GE.1 C C IPROB = 0 NATURAL VIBRATION C IPROB = 1 LINEARIZED BUCKLING ANALYSIS C_______________________________________________________________________ 167 IF(NEIGEN.EQ.0) GO TO 157 DO 710 I = 1,NEQ IF(GSTIF(I,I).EQ.0) GSTIF(I,I) = 1.0 710 CONTINUE C CALL EGNBOU (GSTIF,GMASS,IBDY,NEQ,NEQR,NBDY,NRMAX) C CALL DGVLRG (NEQR,GSTIF,NRMAX,GMASS,NRMAX,ALPHA,EGNVAL) C T11 = REAL(ALPHA1)/EGNVAL(1) C IF(IPROB.EQ.0) THEN C WRITE(IT,840) C DO 276 I = 1,NEQR C FREQ = DSQRT(DABS(EGNVAL(I))) C FREQ1 = TH*FREQ*DSQRT(RHO/E2) C FREQ2 = TH*DSQRT(RHO*T11/E2) C WRITE(IT,845) I,T11,FREQ2 C276 CONTINUE C ELSE C WRITE (IT,750) C WRITE(IT,850) C DO 376 I = 1,NEQR C T22 = REAL(ALPHA1) C WRITE(IT,855) I,EGNVAL(I),T11,T22 C376 CONTINUE C END IF C STOP C____________________________________________________________________ C C.....IMPOSE THE PRIMARY AND SECONDARY B.C'S C____________________________________________________________________ 157 DDLA = DLAMDA(LL)+DLAM DO 713 I = 1,NEQ IF(GSTIF(I,1).EQ.0) GSTIF(I,1) = 1.0 713 CONTINUE IF(NBSF.EQ.0) GO TO 200 DO 170 NF =1,NBSF NB = IBSF(NF) IF(IPOST.EQ.1) THEN FBAR(NB) = FBAR(NB)+DDLA*VBSF(NF) END IF IF(IPOST.NE.1.OR.NITER.EQ.1.OR.IR.EQ.0) THEN GF(NB) = GF(NB)+VBSF(NF) END IF 170 CONTINUE 200 CONTINUE IF(NBDY.EQ.0) GO TO 395 IF(IPOST.EQ.1)THEN CALL BNDY(NRMAX,NCMAX,NEQ,NHBW,GSTIF,FBAR,NBDY,IBDY,VBDY) END IF IF(IPOST.NE.1.OR.NITER.EQ.1.OR.IR.EQ.0) THEN CALL BNDY(NRMAX,NCMAX,NEQ,NHBW,GSTIF,GF,NBDY,IBDY,VBDY) END IF IF(NITER.GT.1.OR.LL.GT.1) GO TO 156 DO 148 IH=1,NBDY 148 VBDY(IH) = 0.0 156 CONTINUE IRES=0 395 CONTINUE IF (IPOST.EQ.0) GO TO 113 C_______________________________________________________________________ C C.... TRACE THE POST-BUCKLING PATH USING THE MODIFIED RIKS ALGORITHM C C IR = 0 use the Full Newton Raphson Method C IR = 1 use Modified Newton Raphson Method C_______________________________________________________________________ ITROPT = 4 CALL RIKS (IT,NRMAX,NCMAX,NEQ,NHBW,IR,ISTF,LL,NITER,ITROPT,ITERP, + ICNVRG,IERR,W,DELU,FBAR,GF,GSTIF,DARC,EPSLON,ERROR) IF(ICNVRG.EQ.0) GO TO 158 GO TO 94 C_______________________________________________________________________ C C.... SOLVE THE BANDED SYMMETRIC EQUATIONS C_______________________________________________________________________ 113 CONTINUE 868 FORMAT(I3,8(E9.2,1X)) CALL SOLVE(NRMAX,NCMAX,NEQ,NHBW,GSTIF,GF,IRES,DET) DO 630 IS=1,NEQ 630 W(IS)= W(IS)+GF(IS) CALL SECOND(T(5)) C______________________________________________________________________ C C.... CHECK FOR CONVERGENCE OF NEWTON-RAPHSON METHOD C_______________________________________________________________________ SUM1 = 0.0 SUM2 = 0.0 DO 99 I=1,NEQ SUM1 = SUM1+GF(I)**2 99 SUM2 = SUM2+W(I)**2 ERROR= DSQRT(SUM1/SUM2) IF(ERROR.GT.EPSLON) GO TO 158 C_______________________________________________________________________ C C..... PRINT OUT DISPLACEMENTS C_______________________________________________________________________ 94 CONTINUE DLAMDA(LL+1) = DLAMDA(LL) + DLAM IF(IPOST.EQ.0) DLAMDA(LL+1) = 1.0 IF(IPOST.EQ.0) DDLA = 1.0 DLOD = 0.0 DO 67 IK = 1,NBSF 67 DLOD = DLOD+VBSF(IK)*DLAMDA(LL+1) IF(NBSF.NE.0) DLOAD = DLOD IF(NBSF.EQ.0) DLOAD = DDLA*PW IF(NBSF.EQ.0.AND.ILOD.EQ.2) DLOAD = 2*DDLA*PW WRITE(6,394) DLOAD WRITE(6,505) WRITE (6,515) IF(ISMEAR.EQ.1) NIF= NPLY+1+NPLS+NPLR DO 637 I = 1,NEQ 637 IF(DABS(W(I)).LT.1E-06) W(I) = 0.0 DO 518 J=1,NNODE J1=(J-1)*NDF+1 DO 528 JJ=1,NIF J2 = J1 + NIF J3 = J1 + 2*NIF IF(JJ.EQ.1) THEN WRITE(IT,538) J,JJ,W(J1),W(J2),W(J3) ELSE WRITE(IT,548) JJ,W(J1),W(J2),W(J3) END IF 528 J1=J1+1 518 CONTINUE IF(IERR.LT.5) ITERP = NITER WRITE(6,98) NITER,ERROR 98 FORMAT(/,5X,'Iteration number =',I2,/,5X,'Measure of Error =', + 1X,F7.5/) C_______________________________________________________________________ C C.... POST PROCESSING OF SOLUTION C C_______________________________________________________________________ C_______________________________________________________________________ C C.... COMPUTE THE DERIVATIVES OF THE SOLUTION C_______________________________________________________________________ 230 CONTINUE IB = (NNODE+1)*NDF/2 IF(NBSF.NE.0) IB = IBSF(1) C IF(IPOST.EQ.0) THEN C WRITE(7,436) LL,DLOAD,W(IB),W(IB+1),W(IB+2) C ELSE C WRITE(7,436) LL,DLOAD,W(IB),W(IB-2),DLAM,DARC,DDLA C 436 FORMAT(3X,I2,6(E10.3,1X)) C END IF IF(LL.EQ.NLOAD.AND.ISTRES.GE.1) THEN C_______________________________________________________________________ C C CALCULATE THE STRAINS AND STRESSES OF THE SHELL C_______________________________________________________________________ ICT = 1 IF (NELEM.GT.5) ICT = 5 DO 240 N = 1,NELEM,ICT L=0 DO 242 I=1,NPE NI=NOD(N,I) LI=(NI-1)*NDF ELXY(I,1)=X(NI) ELXY(I,2)=Y(NI) DO 242 J=1,NDF LI=LI+1 L=L+1 GF(L) = W(LI) 242 CONTINUE C_______________________________________________________________________ 240 CALL STRESS(NPE,NDF,NIF,N,ICODE) END IF C_______________________________________________________________________ C_______________________________________________________________________ C C CALCULATE THE STRAINS AND STRESSES OF THE STIFFENERS C___________________________________________________________ ____________ IF(LL.EQ.NLOAD.AND.ISTRES.EQ.2) THEN WRITE (6,62) IFLAG = 1 RAD = 1E12 DO 775 N = 1,NSTRG,ict L = 0 DO 778 I = 1,NPB NI = NODS(N,I) ELX(I) = X(NI) LI = (NI-1)*NDF DO 778 J = 1,NDF LI = LI+1 L = L+1 GF(L) = W(LI) 778 CONTINUE 775 CALL STRSTF(NPB,NDF,NIF,N,ICODE,RAD,IFLAG,THETS) IFLAG = 2 DO 785 N = 1,NRING,ict IF (ILOCR(N).EQ.0) THEN RAD = RADIUS(1) ELSE RAD = RADIUS(NPLY+1) END IF L = 0 DO 788 I = 1,NPB NI = NODR(N,I) ELX(I) = Y(NI) LI = (NI-1)*NDF DO 788 J = 1,NDF LI = LI+1 L = L+1 GF(L) = W(LI) 788 CONTINUE 785 CALL STRSTF(NPB,NDF,NIF,N,ICODE,RAD,IFLAG,THETR) END IF C_______________________________________________________________________ IF(NLOAD.GT.1.AND.IPOST.EQ.0) THEN READ(IN,*) PX,PW,(VBSF(I),I=1,NBSF) END IF IF(IPOST.EQ.1) THEN DO 451 I = 1,NEQ 451 GUPLS(I) = W(I) END IF 55 CONTINUE C_______________________________________________________________________ C C READ-WRITE FOR RESTART OPERATIONS C C solution W, arc-length DARC, load parameters DLAM and DLAMDA(1), C converged iteration number ITERP are saved for restart analysis. C C ISTART = 0 No Read/Write C ISTART = 1 Write C ISTART = 2 Read C ISTART = 3 Read/Write C_______________________________________________________________________ IF(ISTART.EQ.1.OR.ISTART.EQ.3) THEN REWIND 10 WRITE (10,101) ITERP,DARC,DLAM,DLAMDA(LL) WRITE (10,102) (W(I),I = 1,NEQ) END IF CALL SECOND(T(6)) TT = T(6) - T(1) DO 500 I=1,5 500 T(I)=T(I+1)-T(I) T(6) = TT WRITE (6,1000) (T(i),i=1,6) C_______________________________________________________________________ C C.....F O R M A T S C_______________________________________________________________________ 222 FORMAT(//, 4X,67('_'),//, + 10X,5('C'),4X,4('O'),6X,5('S'),3X,'H',4X,'E',3X, + 6('E'),2X,'L'7X,'L',/, + 9X,('C'),8X,('O'),4X,('O'),4X,'/',8X,'H',4X, + ('H'),3X,'E',7X,'L',7X,'L'/, + 9X,'C',8X,'O',4X,'O',4X,'\',4('S'),'\',3X, + 6('H'),3X,4('E'),4X,'L',7X,'L'/, + 9X,'C',8X,'O',4X,'O',9X,'S',3X,'H',4X, + 'H',3X,'E',7X,'L',7X,'L'/, + 10X,5('C'),4X,4('O'),6X,5('S'),3X,'H',4X,'E',3X, + 6('E'),2X,6('L'),3X,6('L'),//,30x,'Version 1.0',/,29X, + 'December 1992',//,24X,'By Samuel Kinde Kassegne',//) 62 FORMAT(///,5X,18('*'),' Stresses in Stiffeners Follow ',18('*'),/) 300 FORMAT(20A4) 310 FORMAT(/,5X,'Mesh Information :',/,5X,19('_'),/) 315 FORMAT(//,5X,'________ Global Coordinates Follow ______',/, + /,5X,'Node',5X,'X',12X,'Y',/) 394 FORMAT(/,5X,22('_')/,5X,'Load Level =',1PE10.3,/,5X,22('_')/) 325 FORMAT(5X,67('_')) 340 FORMAT(5X,17A4) 345 FORMAT(/,5X,'Num of Elements in the Mesh...... =',I4,/ + ,5X,'Num of Nodes in the Mesh......... =',I4,/ + ,5X,'Num of Nodes per Element ........ =',I4,/ + ,5X,'Num of Layers ................... =',I4,/ + ,5X,'Num of Dof/Node ................. =',I4,/ + ,5X,'Num of Equations ................ =',I4,/ + ,5X,'Half-band Width ................. =',I4,/ + ,5X,'Length of Cylinder ...............=',1PE10.2,/ + ,5X,'Radius of Cylinder ...............=',1PE10.2,/) 348 FORMAT(/,5X,'Code for Order of Integration.....=',I3,// + ,7X,'0 - Full Integration',/ + ,7x,'1 - Reduced (Qi4 & Qi5 terms)',/ + ,7x,'2 - Reduced (Qi3 terms)',/ + ,7x,'3 - Reduced (Qi3, Qi4 & Qi5 terms)',/ + ,7x,'4 - Uniformly Reduced For All Terms',// + ,5X,'Code for Locus of Calc. Stresses =',I3,// + ,7X,'0 - Element Corners ',/ + ,7X,'1 - Gauss Points ',// + ,5X,'Code for Element Stress Computation =',I3,// + ,7X,'0 - No Stress Computations',/ + ,7X,'1 - Stress Computations for Shell Skin',/ + ,7X,'2 - Stress Computations for Shell and Stiffeners',/// + ,5X,'Maximum Number of Iterations ....=',I3,/ + ,5X,'Error Tolerance ...............=',2X,F7.5,/) 353 FORMAT(/,5X,'Code for Type of Analysis........= ',I3,// + ,7X,'1 - Linear Static Analysis',/ + ,7x,'2 - Nonlinear Static Analysis',/ + ,7x,'3 - Postbuckling Analysis',/ + ,7x,'4 - Natural Vibration',/ + ,7x,'5 - Linearized Buckling Analysis',/ + ,7x,'6 - Imperfection Sensitivity Analysis',/) 333 FORMAT(/,5X,'Code for Re-start Operations..... =',I3,// + ,7X,'0 - No Read/Write Activated ',/ + ,7X,'1 - Write Activated ',/ + ,7X,'2 - Read Activated ',/ + ,7X,'3 - Read/Write Activated ',/) 343 FORMAT(/,5X,'Code for Postbuckling Procedure.. =',I3,// + ,7X,'0 - Full Newton-Raphson ',/ + ,7X,'1 - Modified Newton-Raphson ',/) 363 FORMAT(/,5X,'Code for Type of Load............ =',I3,// + ,7X,'0 - Uniformly Distributed Load',/ + ,7x,'1 - Double Sinusoidal Load (both in x-and y-dir.)',/ + ,7x,'2 - Single Sinusoidal Load (only in x-dir)',/ + ,7x,'3 - Single Sinusoidal Load (only in y-dir)',/ + ,7X,'4 - Point Load ',/) 273 FORMAT(/,5X,'Code for Location of Load.......=',I3,// + ,7X,'0 - External pressure',/ + ,7x,'1 - Internal pressure',/ + ,7x,'2 - Uniform pressure (through thickness)',/) 354 FORMAT(/,5X,'Primary Dof Value',/, + 5x,48('_'),/) 357 FORMAT(/,5X,'Secondary Dof Value',/ + 5x,48('_'),/) 360 FORMAT(5X,I4,31X,E12.5) 361 FORMAT(/,5X,'Shell 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,//) 392 FORMAT(/,5X,'Number of Stringers.. = ',I2,/, + 5X,'Number of Rings...... = ',I2,/, + 5X,'Location of Stringer..= ',I2,/, + 5X,'Location of Ring..... = ',I2,/, + 5X,'Mod. of Elasticity(Es)= ',1PE10.3,/, + 5X,'Shear Modulus.....(Gs)= ',1PE10.3,/, + 5X,'Mod.of Elasticity(Es3)= ',1PE10.3,/, + 5X,'Area os Stringer..(As)= ',1PE10.3,/, + 5X,'Mod. of Elasticity(Er)= ',1PE10.3,/, + 5X,'Shear Modulus.....(Gr)= ',1PE10.3,/, + 5X,'Mod.of Elasticity(Er3)= ',1PE10.3,/, + 5X,'Area of Ring .....(Ar)= ',1PE10.3,/) 461 FORMAT(/,5X,'Shell Stiffener Information (Discrete Approach) :', + / ,5X,49('_'),/) 462 FORMAT(/,5X,'Shell Stiffener Information (Smeared Approach) :', + / ,5X,49('_'),/) 463 FORMAT(/,5X,'Shell Stiffener Information (Full Discrete Approach)' + ,/ ,5X,52('_'),/) 492 FORMAT(/,5X,'Spacing of Stringers = ',1PE10.3,/, + 5X,'Eccent. of Stringers = ',1PE10.3,/, + 5X,'Density of Str. Mater.= ',1PE10.3,/, + 5X,'Location of Stringer = ',I2,/, + 5X,'Mod. of Elasticity(ES)= ',1PE10.3,/, + 5X,'Shear Modulus (GS)... = ',1PE10.3,/, + 5X,'X-Area of Stringer... = ',1PE10.3,/, + 5X,'Moment of Inertia(EIS)= ',1PE10.3,/, + 5X,'Polar Mo. Inertia(EJS)= ',1PE10.3,/, + 5x,' ........... ',//, + 5X,'Spacing of Rings..... = ',1PE10.3,/, + 5X,'Eccent. of Rings..... = ',1PE10.3,/, + 5X,'Density of Rng. Mater.= ',1PE10.3,/, + 5X,'Location of Ring..... = ',I2,/, + 5X,'Mod. of Elasticity(ER)= ',1PE10.3,/, + 5X,'Shear Modulus (GR)... = ',1PE10.3,/, + 5X,'X-Area of Ring....... = ',1PE10.3,/, + 5X,'Moment of Inertia(EIR)= ',1PE10.3,/, + 5X,'Polar Mo. Inertia(EJR)= ',1PE10.3,/) 592 FORMAT(/,5X,'Number of Stringers.. = ',I2,/, + 5X,'Number of Rings...... = ',I2,/, + 5X,'Number of Strg.layers = ',I2,/, + 5X,'Number of Ring layers = ',I2,/, + 5X,'Thick. of Strin. Plies= ',F5.3/, + 5X,'Thick. of Ring Plies = ',F5.3/, + 5X,'Mod. of Elasticity(E1)= ',1PE10.3,/, + 5X,'Mod. of Elasticity(E3)= ',1PE10.3,/, + 5X,'Shear Modulus....(G13)= ',1PE10.3) 593 FORMAT(5X,'Location of Stringer..= ',10I2) 594 FORMAT(5X,'Location of Ring..... = ',10I2,/) 371 FORMAT(5X,'Ply Number Ply Angle Ply Thickness', + /,5X,49('_'),/) 373 FORMAT(5X,49('_'),/) 381 FORMAT(5X,I1,12X,F6.2,17X,F6.3) 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, + 'V-j',13X,'W-j',/,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,I3,3X,1PE10.3,3X,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('_'),//) 840 FORMAT(/,15X,28('_'),/,15X, + 'I Eigenvalue Frequency',/,15X,29('_'),/) 845 FORMAT(13X,I3,4X,3(E11.4,2X)) 850 FORMAT(/,5X,58('_'),/,5X,'I',11X,'Pcr',9x,'Egnvec(1)',3x, + 'Egnvec(2)',3x,'Egnvec(3)',/,5X,58('_'),/) 855 FORMAT(3X,I3,8X,E10.3,4X,3(E10.3,2X)) 1000 FORMAT (///' Summary of Solution Time',/,1x,25('_'),// + 41H Form Ply Stiffness ................... = ,F7.2/ + 41H Form Shell Stiffness ................ = ,F7.2/ + 41H Form Stiffener Stiffness ............. = ,F7.2/ + 41H Solve Static Load Case ............... = ,F7.2/ + 41H Compute and Print Stresses and Disp... = ,F7.2// + 41H Total Time Required .................. = ,F7.2,' Seconds'/) STOP END C______________________________________________________________________ C C C... THIS SUBROUTINE CALCULATES THE ELEMENT TANGENT STIFFNESS C... MATRIX OF A LAYERWISE CYLINDRICAL SHELL ELEMENT. C C C______________________________________________________________________ SUBROUTINE STFSTR(PX,PW,NPE,N,RHO,DEP,NDF,NN) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NMX=10,NEX=10,NUX=486,NGX=150,NBX=126) COMMON /SEL/ ST(NGX,NGX),ELP(NGX),ELXY(9,2),ELFC(NGX) COMMON /DOUBLE/ DWDX(20),DWDY(20),TANK(NGX,NGX),DUDX(20),DUDY(20) COMMON /GLOBAL/ GSTIF(NUX,NBX),GF(NUX),W(NUX),EPL(NUX) COMMON /SMAS/ GMASS(NEX,NEX),EGNVAL(NEX) COMMON /ABDF/ A(10,10,6,6), ABAR(10,10,6,6),DBAR(10,10,6,6), + G(10,10) COMMON /DDBF/ D(10,10,10,6,6),DB(10,10,10,6,6),F(10,10,10,10,6,6) COMMON /PAR/ NRED,ILOD,IPROB,IPOST,NIT,LL,IAXIAL,ILATER,ISHEAR COMMON /DWDS/ DVDX(20),DVDY(20),DWKDX(20),DWKDY(20) COMMON /ELAST/ THETA(20),RADIUS(15),ELX(3) COMMON /DEPTH/ Z(15),ZZ(15),THICK(15),SM(NGX,NGX) COMMON /EDWR/ EDWDX(25,20),EDWDY(25,20) COMMON /PAR1/ NEIGEN,IFLAG,NPLY,NPLS,NPLR COMMON /DLDA/ DLAMDA(15),DLAM COMMON /PAR2/ IMPFCT,ISINE,ISMEAR COMMON /SPAND/ SPANX,SPANY DIMENSION WT(4,4),GAUSS(4,4),SF(9),GSF(2,9) DATA GAUSS / 4*0.0D0,-0.57735027D0,0.57735027D0,2*0.0D0, + -0.77459667D0,0.0D0,0.77459667D0,0.0D0, + -0.86113631D0,-0.33998104D0,0.33998104D0, + 0.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 PI= DACOS(-1.0D0) NIF = NPLY+1 NDF = 3*NIF NN = NPE*NDF IF(NPE.EQ.4) THEN NGP = 2 LGP = 2 MGP = 2 IF(NRED.EQ.1.OR.NRED.EQ.3) LGP = 1 IF(NRED.EQ.2.OR.NRED.EQ.3) MGP = 1 ELSE NGP = 3 LGP = 3 MGP = 3 IF(NRED.EQ.1.OR.NRED.EQ.3) LGP = 2 IF(NRED.EQ.2.OR.NRED.EQ.3) MGP = 2 END IF DDLA = DLAMDA(LL)+DLAM IF(IPOST.EQ.0) DDLA = 1.0 DO 5 I = 1,NGX ELP(I) = 0.0 EPL(I) = 0.0 DO 5 J = 1,NGX TANK(I,J) = 0.0 SM(I,J) = 0.0 5 ST(I,J) = 0.0 IF (NRED.EQ.4) THEN NGP = NGP-1 LGP = LGP-1 MGP = MGP-1 END IF KGP = 4 C_______________________________________________________________________ C C.....LOAD VECTOR FOR DOUBLE SINUSOIDAL LOADS FOLLOWS C C ISINE = 0 ; Uniform load C ISINE = 1 ; Double sinusoidal load (both in x and y direction) C ISINE = 2 ; Single sinusoidal load (in x direction direction only) C ISINE = 3 ; Single sinusoidal load (in y direction direction only) C C.... ILOD = 2 ; Both internal and external sinusoidal load C.... ILOD = 1 ; Internal load only. C.... ILOD = 0 ; External load only. C C_______________________________________________________________________ IF(ISINE.GT.0) THEN DO 22 NI = 1,KGP DO 22 NJ = 1,KGP XI = GAUSS(NI,KGP) ETA= GAUSS(NJ,KGP) CALL SHAPE2(NPE,XI,ETA,SF,GSF,DET,ELXY) CNST = DET*WT(NI,KGP)*WT(NJ,KGP) XX = 0.0 YY = 0.0 DO 28 I = 1,NPE XX = XX + SF(I)*ELXY(I,1) 28 YY = YY + SF(I)*ELXY(I,2) DO 32 I = 1,NPE IF(ISINE.EQ.1) THEN PWT = PW*DSIN(PI*XX/SPANX)*DSIN(PI*YY/SPANY)*SF(I)*CNST ELSE IF(ISINE.EQ.2) THEN PWT = PW*DSIN(PI*XX/SPANX)*SF(I)*CNST ELSE IF(ISINE.EQ.3) THEN PWT = PW*DSIN(PI*YY/SPANY)*SF(I)*CNST END IF L1 = (I-1)*NDF + 1 IF(ILOD.EQ.2) GO TO 114 IF(ILOD.EQ.0) THEN I2 = L1 + 3*NIF-1 ELSE I2 = L1 + 2*NIF END IF ELP(I2) = ELP(I2) + PWT GO TO 139 114 DO 138 I1 = L1, L1+NIF-1 I2 = I1 + 2*NIF 138 ELP(I2) = ELP(I2) + PWT 139 CONTINUE 32 CONTINUE 22 CONTINUE END IF C_______________________________________________________________________ C C......FULL INTEGRATION OF TANGENT STIFFNESS MATRIX FOLLOWS C C_______________________________________________________________________ DO 100 NI =1,NGP DO 100 NJ =1,NGP XI = GAUSS(NI,NGP) ETA= GAUSS(NJ,NGP) CALL SHAPE2(NPE,XI,ETA,SF,GSF,DET,ELXY) CNST = DET*WT(NI,NGP)*WT(NJ,NGP) C C.... LOAD VECTOR FOLLOWS C IF(ISINE.EQ.0) THEN DO 110 I=1,NPE L1 = (I-1)*NDF +1 C C.... ILOD = 2 ; Both internal and external pressure. C.... ILOD = 1 ; Internal pressure only. C.... ILOD = 0 ; External pressure only. C IF(ILOD.EQ.2) GO TO 115 IF(ILOD.EQ.0) THEN I1 = L1+NIF-1 I2 = L1 + 3*NIF-1 ELSE I1 = L1 I2 = L1 + 2*NIF END IF ELP(I1) = ELP(I1) + DDLA*PX*SF(I)*CNST ELP(I2) = ELP(I2) + DDLA*PW*SF(I)*CNST EPL(I1) = EPL(I1) + PX*SF(I)*CNST EPL(I2) = EPL(I2) + PW*SF(I)*CNST GO TO 110 115 DO 118 I1=L1,L1+NIF-1 I2=I1+2*NIF ELP(I1) = ELP(I1) + DDLA*PX*SF(I)*CNST ELP(I2) = ELP(I2) + DDLA*PW*SF(I)*CNST EPL(I1) = EPL(I1) + PX*SF(I)*CNST 118 EPL(I2) = EPL(I2) + PW*SF(I)*CNST 110 CONTINUE END IF DO 15 I = 1,20 DUDX(I) = 0.0 DUDY(I) = 0.0 DVDX(I) = 0.0 DVDY(I) = 0.0 DWDX(I) = 0.0 DWDY(I) = 0.0 DWKDX(I) = 0.0 15 DWKDY(I) = 0.0 DO 130 I=1,NPE L1 =(I-1)*NDF+2*NIF+1 L2 =(I-1)*NDF+1 L3 =(I-1)*NDF+NIF+1 L4 =(I-1)*NDF+3*NIF+1 DO 140 J= 1,NIF DUDX(J) = DUDX(J)+GSF(1,I)*ELFC(L2) DUDY(J) = DUDY(J)+GSF(2,I)*ELFC(L2) DVDX(J) = DVDX(J)+GSF(1,I)*ELFC(L3) DVDY(J) = DVDY(J)+GSF(2,I)*ELFC(L3) DWDX(J) = DWDX(J)+GSF(1,I)*ELFC(L1) DWDY(J) = DWDY(J)+GSF(2,I)*ELFC(L1) EDWDX(N,J) = DWDX(J) EDWDY(N,J) = DWDY(J) IF(IMPFCT.EQ.1) THEN DWKDX(J) = DWKDX(J)+GSF(1,I)*ELFC(L4) DWKDY(J) = DWKDY(J)+GSF(2,I)*ELFC(L4) END IF L1 = L1+1 L2 = L2+1 L3 = L3+1 L4 = L4+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 IF(IPROB.EQ.0) THEN CNST1 = RHO*SF(I)*SF(J)*CNST CNST2 = CNST1 ELSE CNST2=(IAXIAL*GSF(1,I)*GSF(1,J) + + ISHEAR*(GSF(1,I)*GSF(2,J)+GSF(2,I)*GSF(1,J)) + + ILATER*GSF(2,I)*GSF(2,J))*CNST/DEP CNST1 = 0.0 END IF DO 55 IN = JJ,JN NE = NE+1 R = RADIUS(NE) CON11 = 0.0 CON22 = 0.0 CON22I = 0.0 CON33 = 0.0 CON44 = 0.0 CON55 = 0.0 CONAA = 0.0 CONBB = 0.0 CON66 = 0.0 CON66I = 0.0 CON555 = 0.0 CON666 = 0.0 CON311 = 0.0 CON312 = 0.0 CON321 = 0.0 CON322 = 0.0 CON77 = 0.0 CON88 = 0.0 CON99 = 0.0 CON97 = 0.0 CON97I = 0.0 CON14 = 0.0 CON141 = 0.0 CON14I = 0.0 CON24 = 0.0 CON242 = 0.0 CON34 = 0.0 CON341 = 0.0 CON34I = 0.0 CON344 = 0.0 CON355 = 0.0 CON34A = 0.0 CON34B = 0.0 CON34C = 0.0 CON1D = 0.0 CON2D = 0.0 CON22D = 0.0 CON3D = 0.0 CON44D = 0.0 CON55T = 0.0 CON55S = 0.0 CON66T = 0.0 CON66S = 0.0 CON55M = 0.0 CON66M = 0.0 CON66J = 0.0 C_______________________________________________________________________ C C.....CALCULATE [K11] (Tangent and Direct) and [S11] C C_______________________________________________________________________ ST(I1,IN) = ST(I1,IN)+(A(NE,NE,1,1)*GSF(1,I)*GSF(1,J) + + A(NE,NE,1,6)*GSF(1,I)*GSF(2,J) + + A(NE,NE,1,6)*GSF(1,J)*GSF(2,I) + + A(NE,NE,6,6)*GSF(2,I)*GSF(2,J))*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(1,I)*GSF(1,J) + + A(NE,NE+1,1,6)*GSF(1,I)*GSF(2,J) + + A(NE,NE+1,1,6)*GSF(1,J)*GSF(2,I) + + A(NE,NE+1,6,6)*GSF(2,I)*GSF(2,J))*CNST ST(I1+1,IN) = ST(I1,IN+1) TANK(I1,IN+1) = ST(I1,IN+1) TANK(I1+1,IN) = ST(I1,IN+1) 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 ST(I1,IN+NIF) = ST(I1,IN+NIF) + +(A(NE,NE,1,2)*GSF(1,I)*GSF(2,J) + + A(NE,NE,1,6)*GSF(1,J)*GSF(1,I) + + A(NE,NE,2,6)*GSF(2,I)*GSF(2,J) + + A(NE,NE,6,6)*GSF(2,I)*GSF(1,J))*CNST TANK(I1,IN+NIF) = ST(I1,IN+NIF) IF(IN.LT.JN) THEN ST(I1,IN+NIF+1) = ST(I1,IN+NIF+1) + +(A(NE,NE+1,1,2)*GSF(1,I)*GSF(2,J) + + A(NE,NE+1,1,6)*GSF(1,J)*GSF(1,I) + + A(NE,NE+1,2,6)*GSF(2,I)*GSF(2,J) + + A(NE,NE+1,6,6)*GSF(2,I)*GSF(1,J))*CNST ST(I1+1,IN+NIF) = ST(I1,IN+NIF+1) TANK(I1,IN+NIF+1) = ST(I1,IN+NIF+1) TANK(I1+1,IN+NIF) = ST(I1,IN+NIF+1) END IF C_______________________________________________________________________ C C.....CALCULATE [K21] (Tangent and Direct) C C_______________________________________________________________________ C ST(I1+NIF,IN) = ST(I1+NIF,IN) + +(A(NE,NE,1,2)*GSF(2,I)*GSF(1,J) + + A(NE,NE,1,6)*GSF(1,I)*GSF(1,J) + + A(NE,NE,2,6)*GSF(2,I)*GSF(2,J) + + A(NE,NE,6,6)*GSF(1,I)*GSF(2,J))*CNST TANK(I1+NIF,IN) = ST(I1+NIF,IN) IF(IN.LT.JN) THEN ST(I1+NIF,IN+1) = ST(I1+NIF,IN+1) + +(A(NE,NE+1,1,2)*GSF(2,I)*GSF(1,J) + + A(NE,NE+1,1,6)*GSF(1,I)*GSF(1,J) + + A(NE,NE+1,2,6)*GSF(2,I)*GSF(2,J) + + A(NE,NE+1,6,6)*GSF(1,I)*GSF(2,J))*CNST ST(I1+NIF+1,IN) = ST(I1+NIF,IN+1) TANK(I1+NIF,IN+1) = ST(I1+NIF,IN+1) TANK(I1+NIF+1,IN) = ST(I1+NIF,IN+1) END IF C_______________________________________________________________________ C C.....CALCULATE [K22] (Tangent and Direct) and [S22] C C_______________________________________________________________________ C ST(I1+NIF,IN+NIF) = ST(I1+NIF,IN+NIF) + +(A(NE,NE,2,2)*GSF(2,I)*GSF(2,J) + + A(NE,NE,2,6)*GSF(1,J)*GSF(2,I) + + A(NE,NE,2,6)*GSF(1,I)*GSF(2,J) + + A(NE,NE,6,6)*GSF(1,J)*GSF(1,I))*CNST TANK(I1+NIF,IN+NIF) = ST(I1+NIF,IN+NIF) SM(I1+NIF,IN+NIF) = SM(I1+NIF,IN+NIF)+G(NE,NE)*CNST1 IF(IN.LT.JN) THEN ST(I1+NIF,IN+NIF+1) = ST(I1+NIF,IN+NIF+1) + +(A(NE,NE+1,2,2)*GSF(2,I)*GSF(2,J) + +A(NE,NE+1,2,6)*GSF(1,J)*GSF(2,I) + +A(NE,NE+1,2,6)*GSF(1,I)*GSF(2,J) + +A(NE,NE+1,6,6)*GSF(1,J)*GSF(1,I))*CNST ST(I1+NIF+1,IN+NIF) = ST(I1+NIF,IN+NIF+1) TANK(I1+NIF,IN+NIF+1) = ST(I1+NIF,IN+NIF+1) TANK(I1+NIF+1,IN+NIF) = ST(I1+NIF,IN+NIF+1) SM(I1+NIF,IN+NIF+1) = SM(I1+NIF,IN+NIF+1)+G(NE,NE+1)*CNST1 SM(I1+NIF+1,IN+NIF) = SM(I1+NIF,IN+NIF+1) END IF C C_______________________________________________________________________ C C CALCULATE [K13] (Tangent and Direct) C C_______________________________________________________________________ C C C ....LINEAR TERMS FOLLOW C IP = IN+2*NIF CON1 = A(NE,NE,1,2)*GSF(1,I)*SF(J)/R + +A(NE,NE,2,6)*GSF(2,I)*SF(J)/R IF(IN.LT.JN) THEN CON2 = A(NE,NE+1,1,2)*GSF(1,I)*SF(J)/R + + A(NE,NE+1,2,6)*GSF(2,I)*SF(J)/R CON2I= A(NE+1,NE,1,2)*GSF(1,I)*SF(J)/R + + A(NE+1,NE,2,6)*GSF(2,I)*SF(J)/R END IF C C.... NON-LINEAR TERMS OF K13 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(1,I)*GSF(1,J)*(DWDX(ME)+DWKDX(ME)) + +2*D(NE,NE,ME,1,6)*GSF(1,I)*GSF(1,J)*(DWDY(ME)+DWKDY(ME)) + +D(NE,NE,ME,1,2)*GSF(1,I)*GSF(2,J)*(DWDY(ME)+DWKDY(ME)) + +D(NE,NE,ME,1,6)*GSF(2,I)*GSF(1,J)*(DWDX(ME)+DWKDX(ME)) + +D(NE,NE,ME,2,6)*GSF(2,I)*GSF(2,J)*(DWDY(ME)+DWKDY(ME)) + +2*D(NE,NE,ME,6,6)*GSF(2,I)*GSF(1,J)*(DWDY(ME)+DWKDY(ME)) CON1D = CON1D + + D(NE,NE,ME,1,1)*GSF(1,I)*GSF(1,J)*(DWDX(ME)+2*DWKDX(ME))/2 + +D(NE,NE,ME,1,6)*GSF(1,I)*GSF(1,J)*(DWDY(ME)+DWKDY(ME)) + +D(NE,NE,ME,1,2)*GSF(1,I)*GSF(2,J)*(DWDY(ME)+2*DWKDY(ME))/2 + +D(NE,NE,ME,1,6)*GSF(2,I)*GSF(1,J)*(DWDX(ME)+2*DWKDX(ME))/2 + +D(NE,NE,ME,2,6)*GSF(2,I)*GSF(2,J)*(DWDY(ME)+2*DWKDY(ME))/2 + +D(NE,NE,ME,6,6)*GSF(2,I)*GSF(1,J)*(DWDY(ME)+DWKDY(ME)) ME = ME+1 150 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 IX = NE DO 160 IC=1,2 CON22 = CON22 + + D(NE,NE+1,IZ,1,1)*GSF(1,I)*GSF(1,J)*(DWDX(IZ)+DWKDX(IZ)) + +2*D(NE,NE+1,IZ,1,6)*GSF(1,I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) + +D(NE,NE+1,IZ,1,2)*GSF(1,I)*GSF(2,J)*(DWDY(IZ)+DWKDY(IZ)) + +D(NE,NE+1,IZ,1,6)*GSF(2,I)*GSF(1,J)*(DWDX(IZ)+DWKDX(IZ)) + +D(NE,NE+1,IZ,2,6)*GSF(2,I)*GSF(2,J)*(DWDY(IZ)+DWKDY(IZ)) + +2*D(NE,NE+1,IZ,6,6)*GSF(2,I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) CON2D = CON2D + + D(NE,NE+1,IZ,1,1)*GSF(1,I)*GSF(1,J)*(DWDX(IZ)+2*DWKDX(IZ))/2 + +D(NE,NE+1,IZ,1,6)*GSF(1,I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) + +D(NE,NE+1,IZ,1,2)*GSF(1,I)*GSF(2,J)*(DWDY(IZ)+2*DWKDY(IZ))/2 + +D(NE,NE+1,IZ,1,6)*GSF(2,I)*GSF(1,J)*(DWDX(IZ)+2*DWKDX(IZ))/2 + +D(NE,NE+1,IZ,2,6)*GSF(2,I)*GSF(2,J)*(DWDY(IZ)+2*DWKDY(IZ))/2 + +D(NE,NE+1,IZ,6,6)*GSF(2,I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) CON22I= CON22I + + D(NE+1,NE,IZ,1,1)*GSF(1,I)*GSF(1,J)*(DWDX(IZ)+DWKDX(IZ)) + +2*D(NE+1,NE,IZ,1,6)*GSF(1,I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) + +D(NE+1,NE,IZ,1,2)*GSF(1,I)*GSF(2,J)*(DWDY(IZ)+DWKDY(IZ)) + +D(NE+1,NE,IZ,1,6)*GSF(2,I)*GSF(1,J)*(DWDX(IZ)+DWKDX(IZ)) + +D(NE+1,NE,IZ,2,6)*GSF(2,I)*GSF(2,J)*(DWDY(IZ)+DWKDY(IZ)) + +2*D(NE+1,NE,IZ,6,6)*GSF(2,I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) CON22D = CON22D + + D(NE+1,NE,IZ,1,1)*GSF(1,I)*GSF(1,J)*(DWDX(IZ)+2*DWKDX(IZ))/2 + +D(NE+1,NE,IZ,1,6)*GSF(1,I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) + +D(NE+1,NE,IZ,1,2)*GSF(1,I)*GSF(2,J)*(DWDY(IZ)+2*DWKDY(IZ))/2 + +D(NE+1,NE,IZ,1,6)*GSF(2,I)*GSF(1,J)*(DWDX(IZ)+2*DWKDX(IZ))/2 + +D(NE+1,NE,IZ,2,6)*GSF(2,I)*GSF(2,J)*(DWDY(IZ)+2*DWKDY(IZ))/2 + +D(NE+1,NE,IZ,6,6)*GSF(2,I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) IZ = IZ+1 160 CONTINUE END IF ST(I1,IP) = ST(I1,IP)+(CON1+CON1D)*CNST TANK(I1,IP) = TANK(I1,IP)+(CON1+CON11)*CNST IF(IN.LT.JN) THEN ST(I1,IP+1) = ST(I1,IP+1)+(CON2+CON2D)*CNST ST(I1+1,IP) = ST(I1+1,IP)+(CON2I+CON22D)*CNST TANK(I1,IP+1)= TANK(I1,IP+1)+(CON2+CON22)*CNST TANK(I1+1,IP)= TANK(I1+1,IP)+(CON2I+CON22I)*CNST END IF C C_______________________________________________________________________ C C CALCULATE [K31] (Tangent and Direct) C C_______________________________________________________________________ C C C ....LINEAR TERMS FOLLOW C IP = IN IQ = I1+2*NIF CON1 = A(NE,NE,1,2)*GSF(1,J)*SF(I)/R + + A(NE,NE,2,6)*GSF(2,J)*SF(I)/R IF(IN.LT.JN) THEN CON2 = A(NE,NE+1,1,2)*GSF(1,J)*SF(I)/R + + A(NE,NE+1,2,6)*GSF(2,J)*SF(I)/R CON2I= A(NE+1,NE,1,2)*GSF(1,J)*SF(I)/R + +A(NE+1,NE,2,6)*GSF(2,J)*SF(I)/R END IF C C.... NON-LINEAR TERMS OF K31 FOLLOW C ME = NE-1 DO 350 IC=1,3 IF(ME.EQ.0) ME = NE CON99 =CON99 + + D(NE,NE,ME,1,1)*GSF(1,J)*GSF(1,I)*(DWDX(ME)+DWKDX(ME)) + +D(NE,NE,ME,1,2)*GSF(1,J)*GSF(2,I)*(DWDY(ME)+DWKDY(ME)) + +D(NE,NE,ME,1,6)*GSF(2,J)*GSF(1,I)*(DWDX(ME)+DWKDX(ME)) + +D(NE,NE,ME,2,6)*GSF(2,J)*GSF(2,I)*(DWDY(ME)+DWKDY(ME)) + +D(NE,NE,ME,1,6)*GSF(1,J)*(GSF(2,I)*(DWDX(ME)+DWKDX(ME)) + +GSF(1,I)*(DWDY(ME)+DWKDY(ME))) + +D(NE,NE,ME,6,6)*GSF(2,J)*(GSF(2,I)*(DWDX(ME)+DWKDX(ME))+ + GSF(1,I)*(DWDY(ME)+DWKDY(ME))) 344 ME = ME+1 350 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 360 IC=1,2 CON97 = CON97+ + D(NE,NE+1,IZ,1,1)*GSF(1,J)*GSF(1,I)*(DWDX(IZ)+DWKDX(IZ)) + + D(NE,NE+1,IZ,1,2)*GSF(1,J)*GSF(2,I)*(DWDY(IZ)+DWKDY(IZ)) + + D(NE,NE+1,IZ,1,6)*GSF(2,J)*GSF(1,I)*(DWDX(IZ)+DWKDX(IZ)) + + D(NE,NE+1,IZ,2,6)*GSF(2,J)*GSF(2,I)*(DWDY(IZ)+DWKDY(IZ)) + + D(NE,NE+1,IZ,1,6)*GSF(1,J)*(GSF(2,I)*(DWDX(IZ)+DWKDX(IZ))+ + GSF(1,I)*(DWDY(IZ)+DWKDY(IZ))) + + D(NE,NE+1,IZ,6,6)*GSF(2,J)*(GSF(2,I)*(DWDX(IZ)+DWKDX(IZ))+ + GSF(1,I)*(DWDY(IZ)+DWKDY(IZ))) CON97I= CON97I+ + D(NE+1,NE,IZ,1,1)*GSF(1,J)*GSF(1,I)*(DWDX(IZ)+DWKDX(IZ))+ + D(NE+1,NE,IZ,1,2)*GSF(1,J)*GSF(2,I)*(DWDY(IZ)+DWKDY(IZ))+ + D(NE+1,NE,IZ,1,6)*GSF(2,J)*GSF(1,I)*(DWDX(IZ)+DWKDX(IZ))+ + D(NE+1,NE,IZ,2,6)*GSF(2,J)*GSF(2,I)*(DWDY(IZ)+DWKDY(IZ))+ + D(NE+1,NE,IZ,1,6)*GSF(1,J)*(GSF(2,I)*(DWDX(IZ)+DWKDX(IZ))+ + GSF(1,I)*(DWDY(IZ)+DWKDY(IZ)))+ + D(NE+1,NE,IZ,6,6)*GSF(2,J)*(GSF(2,I)*(DWDX(IZ)+DWKDX(IZ))+ + GSF(1,I)*(DWDY(IZ)+DWKDY(IZ))) IZ = IZ+1 360 CONTINUE END IF ST(IQ,IP) = ST(IQ,IP) + (CON1+CON99)*CNST TANK(IQ,IP) = TANK(IQ,IP) + (CON1+CON99)*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1) + (CON2+CON97)*CNST ST(IQ+1,IP) = ST(IQ+1,IP) + (CON2I+CON97I)*CNST TANK(IQ,IP+1) = TANK(IQ,IP+1)+ (CON2+CON97)*CNST TANK(IQ+1,IP) = TANK(IQ+1,IP)+ (CON2I+CON97I)*CNST END IF C C_______________________________________________________________________ C C CALCULATE [K23] (Tangent and Direct) C C_______________________________________________________________________ C C C ....LINEAR TERMS FOLLOW C IP = IN+2*NIF IQ = I1+NIF CON3 = (A(NE,NE,2,2)*GSF(2,I)*SF(J)/R)+ + (A(NE,NE,2,6)*GSF(1,I)*SF(J)/R) IF(IN.LT.JN) THEN CON4 = (A(NE,NE+1,2,2)*GSF(2,I)*SF(J)/R) + +(A(NE,NE+1,2,6)*GSF(1,I)*SF(J)/R) CON4I= (A(NE+1,NE,2,2)*GSF(2,I)*SF(J)/R)+ + (A(NE+1,NE,2,6)*GSF(1,I)*SF(J)/R) END IF C C.... NON-LINEAR TERMS OF K23 FOLLOW C ME = NE-1 DO 170 IC=1,3 IF(ME.EQ.0) ME = NE CON33 = CON33+ + D(NE,NE,ME,1,2)*GSF(2,I)*GSF(1,J)*(DWDX(ME)+DWKDX(ME)) + +2*D(NE,NE,ME,2,6)*GSF(2,I)*GSF(1,J)*(DWDY(ME)+DWKDY(ME)) + +D(NE,NE,ME,2,2)*GSF(2,I)*GSF(2,J)*(DWDY(ME)+DWKDY(ME)) + +D(NE,NE,ME,1,6)*GSF(1,I)*GSF(1,J)*(DWDX(ME)+DWKDX(ME)) + +D(NE,NE,ME,2,6)*GSF(1,I)*GSF(2,J)*(DWDY(ME)+DWKDY(ME)) + +2*D(NE,NE,ME,6,6)*GSF(1,I)*GSF(1,J)*(DWDY(ME)+DWKDY(ME)) CON3D = CON3D+ + D(NE,NE,ME,1,2)*GSF(2,I)*GSF(1,J)*(DWDX(ME)+2*DWKDX(ME))/2 + +D(NE,NE,ME,2,6)*GSF(2,I)*GSF(1,J)*(DWDY(ME)+DWKDY(ME)) + +D(NE,NE,ME,2,2)*GSF(2,I)*GSF(2,J)*(DWDY(ME)+2*DWKDY(ME))/2 + +D(NE,NE,ME,1,6)*GSF(1,I)*GSF(1,J)*(DWDX(ME)+2*DWKDX(ME))/2 + +D(NE,NE,ME,2,6)*GSF(1,I)*GSF(2,J)*(DWDY(ME)+2*DWKDY(ME))/2 + +D(NE,NE,ME,6,6)*GSF(1,I)*GSF(1,J)*(DWDY(ME)+DWKDY(ME)) 165 ME = ME+1 170 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 180 IC=1,2 CON44 = CON44 + + D(NE,NE+1,IZ,1,2)*GSF(2,I)*GSF(1,J)*(DWDX(IZ)+DWKDX(IZ)) + +2*D(NE,NE+1,IZ,2,6)*GSF(2,I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) + +D(NE,NE+1,IZ,2,2)*GSF(2,I)*GSF(2,J)*(DWDY(IZ)+DWKDY(IZ)) + +D(NE,NE+1,IZ,1,6)*GSF(1,I)*GSF(1,J)*(DWDX(IZ)+DWKDX(IZ)) + +D(NE,NE+1,IZ,2,6)*GSF(1,I)*GSF(2,J)*(DWDY(IZ)+DWKDY(IZ)) + +2*D(NE,NE+1,IZ,6,6)*GSF(1,I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) CON44D= CON44D+ + D(NE,NE+1,IZ,1,2)*GSF(2,I)*GSF(1,J)*(DWDX(IZ)+DWKDX(IZ))/2 + +D(NE,NE+1,IZ,2,6)*GSF(2,I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) + +D(NE,NE+1,IZ,2,2)*GSF(2,I)*GSF(2,J)*(DWDY(IZ)+DWKDY(IZ))/2 + +D(NE,NE+1,IZ,1,6)*GSF(1,I)*GSF(1,J)*(DWDX(IZ)+DWKDX(IZ))/2 + +D(NE,NE+1,IZ,2,6)*GSF(1,I)*GSF(2,J)*(DWDY(IZ)+DWKDY(IZ))/2 + +D(NE,NE+1,IZ,6,6)*GSF(1,I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) IZ = IZ+1 180 CONTINUE END IF ST(IQ,IP) = ST(IQ,IP) + (CON3+CON3D)*CNST TANK(IQ,IP) = TANK(IQ,IP) + (CON3+CON33)*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1) + (CON4+CON44D)*CNST ST(IQ+1,IP) = ST(IQ+1,IP) + (CON4I+CON44D)*CNST TANK(IQ,IP+1) = TANK(IQ,IP+1) + (CON4+CON44)*CNST TANK(IQ+1,IP) = TANK(IQ+1,IP) + (CON4I+CON44)*CNST END IF C C_______________________________________________________________________ C C CALCULATE [K32] (Tangent and Direct) C C_______________________________________________________________________ C C C ....LINEAR TERMS FOLLOW C IQ = I1+2*NIF IP = IN+NIF CON3 = A(NE,NE,2,2)*GSF(2,J)*SF(I)/R + + A(NE,NE,2,6)*GSF(1,J)*SF(I)/R IF(IN.LT.JN) THEN CON4 = A(NE,NE+1,2,2)*GSF(2,J)*SF(I)/R + + A(NE,NE+1,2,6)*GSF(1,J)*SF(I)/R CON4I= A(NE+1,NE,2,2)*GSF(2,J)*SF(I)/R + +A(NE+1,NE,2,6)*GSF(1,J)*SF(I)/R END IF C C.... NON-LINEAR TERMS OF K32 FOLLOW C ME = NE-1 DO 370 IC=1,3 IF(ME.EQ.0) ME = NE CON77 = CON77 + + D(NE,NE,ME,1,2)*GSF(2,J)*GSF(1,I)*(DWDX(ME)+DWKDX(ME)) + +D(NE,NE,ME,2,2)*GSF(2,J)*GSF(2,I)*(DWDY(ME)+DWKDY(ME)) + +D(NE,NE,ME,1,6)*GSF(1,J)*GSF(1,I)*(DWDX(ME)+DWKDX(ME)) + +D(NE,NE,ME,2,6)*GSF(1,J)*GSF(2,I)*(DWDY(ME)+DWKDY(ME)) + +D(NE,NE,ME,2,6)*GSF(2,J)*(GSF(2,I)*(DWDX(ME)+DWKDX(ME)) + +GSF(1,I)*(DWDY(ME)+DWKDY(ME))) + +D(NE,NE,ME,6,6)*GSF(1,J)*(GSF(2,I)*(DWDX(ME)+DWKDX(ME)) + +GSF(1,I)*(DWDY(ME)+DWKDY(ME))) 365 ME = ME+1 370 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 380 IC=1,2 CON88 = CON88 + + D(NE,NE+1,IZ,1,2)*GSF(2,J)*GSF(1,I)*(DWDX(IZ)+DWKDX(IZ)) + +D(NE,NE+1,IZ,2,2)*GSF(2,J)*GSF(2,I)*(DWDY(IZ)+DWKDY(IZ)) + +D(NE,NE+1,IZ,1,6)*GSF(1,J)*GSF(1,I)*(DWDX(IZ)+DWKDX(IZ)) + +D(NE,NE+1,IZ,2,6)*GSF(1,J)*GSF(2,I)*(DWDY(IZ)+DWKDY(IZ)) + +D(NE,NE+1,IZ,2,6)*GSF(2,J)*(GSF(2,I)*(DWDX(IZ)+DWKDX(IZ)) + +GSF(1,I)*(DWDY(IZ)+DWKDY(IZ))) + +D(NE,NE+1,IZ,6,6)*GSF(1,J)*(GSF(2,I)*(DWDX(IZ)+DWKDX(IZ)) + +GSF(1,I)*(DWDY(IZ)+DWKDY(IZ))) IZ = IZ+1 380 CONTINUE END IF ST(IQ,IP) = ST(IQ,IP) + (CON3+CON77)*CNST TANK(IQ,IP) = TANK(IQ,IP) + (CON3+CON77)*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1) + (CON4+CON88)*CNST ST(IQ+1,IP) = ST(IQ+1,IP) + (CON4I+CON88)*CNST TANK(IQ,IP+1) = TANK(IQ,IP+1) + (CON4+CON88)*CNST TANK(IQ+1,IP) = TANK(IQ+1,IP) + (CON4I+CON88)*CNST END IF C C_______________________________________________________________________ C C CALCULATE [K33] (Tangent and Direct) and [S33] C C_______________________________________________________________________ C C C ....LINEAR TERMS FOLLOW C IQ = I1+2*NIF IP = IN+2*NIF CON5 = A(NE,NE,2,2)*SF(J)*SF(I)/R/R + + IPROB*PW*GSF(2,I)*GSF(2,J)*G(NE,NE)/DEP SM(IQ,IP) = SM(IQ,IP) + G(NE,NE)*CNST2 IF(IN.LT.JN) THEN CON6 = A(NE,NE+1,2,2)*SF(J)*SF(I)/R/R + +IPROB*PW*GSF(2,I)*GSF(2,J)*G(NE,NE+1)/DEP SM(IQ,IP+1) = SM(IQ,IP+1) + G(NE,NE+1)*CNST2 SM(IQ+1,IP) = SM(IQ,IP+1) END IF C C.... SINGLE NON-LINEAR TERMS OF K33 FOLLOW C ME = NE-1 DO 190 IC=1,3 IF(ME.EQ.0) ME = NE CON55 = CON55 + + D(NE,NE,ME,1,2)*GSF(1,J)*SF(I)*(DWDX(ME)+2*DWKDX(ME))/2/R + +D(NE,NE,ME,2,6)*SF(I)*GSF(1,J)*(DWDY(ME)+DWKDY(ME))/R + +D(NE,NE,ME,2,2)*GSF(2,J)*SF(I)*(DWDY(ME)+2*DWKDY(ME))/2/R + +DB(NE,ME,NE,1,3)*SF(I)*GSF(1,J)*(DWDX(ME)+2*DWKDX(ME))/2 + +DB(NE,ME,NE,2,3)*GSF(2,J)*SF(I)*(DWDY(ME)+2*DWKDY(ME))/2 + +DB(NE,ME,NE,3,6)*SF(I)*GSF(1,J)*(DWDY(ME)+DWKDY(ME)) + +D(NE,NE,ME,1,2)*GSF(1,I)*SF(J)*(DWDX(ME)+DWKDX(NE))/R + +D(NE,NE,ME,2,2)*GSF(2,I)*SF(J)*(DWDY(ME)+DWKDY(NE))/R + +DB(NE,NE,ME,1,3)*GSF(1,I)*SF(J)*(DWDX(ME)+DWKDX(NE)) + +DB(NE,NE,ME,2,3)*GSF(2,I)*SF(J)*(DWDY(ME)+DWKDY(NE)) + +DB(NE,NE,ME,3,6)*SF(J)*(GSF(1,I)*(DWDY(ME)+DWKDY(ME)) + +GSF(2,I)*(DWDX(ME)+DWKDX(ME))) + + D(NE,NE,ME,2,6)*SF(J)*(GSF(1,I)*(DWDY(ME)+DWKDY(ME)) + +GSF(2,I)*(DWDX(ME)+DWKDX(ME)))/R ME = ME+1 190 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 200 IC=1,2 CON66 = CON66 + + D(NE,NE+1,IZ,1,2)*GSF(1,J)*SF(I)*(DWDX(IZ)+2*DWKDX(IZ))/2/R + +D(NE,NE+1,IZ,2,6)*SF(I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ))/R + +D(NE,NE+1,IZ,2,2)*GSF(2,J)*SF(I)*(DWDY(IZ)+2*DWKDY(IZ))/2/R + +DB(NE+1,IZ,NE,1,3)*SF(I)*GSF(1,J)*(DWDX(IZ)+2*DWKDX(IZ))/2 + +DB(NE+1,IZ,NE,2,3)*GSF(2,J)*SF(I)*(DWDY(IZ)+2*DWKDY(IZ))/2 + +DB(NE+1,IZ,NE,3,6)*SF(I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) + +D(NE,NE+1,IZ,1,2)*GSF(1,I)*SF(J)*(DWDX(IZ)+DWKDX(IZ))/R + +D(NE,NE+1,IZ,2,2)*GSF(2,I)*SF(J)*(DWDY(IZ)+DWKDY(IZ))/R + +DB(NE,NE+1,IZ,1,3)*GSF(1,I)*SF(J)*(DWDX(IZ)+DWKDX(IZ)) + +DB(NE,NE+1,IZ,2,3)*GSF(2,I)*SF(J)*(DWDY(IZ)+DWKDY(IZ)) + +DB(NE,NE+1,IZ,3,6)*SF(J)*(GSF(1,I)*(DWDY(IZ)+DWKDY(IZ)) + +GSF(2,I)*(DWDX(IZ)+DWKDX(IZ))) + +D(NE,NE+1,IZ,2,6)*SF(J)*(GSF(1,I)*(DWDY(IZ)+DWKDY(IZ)) + +GSF(2,I)*(DWDX(IZ)+DWKDX(IZ)))/R CON66I = CON66I + + D(NE+1,NE,IZ,1,2)*GSF(1,J)*SF(I)*(DWDX(IZ)+2*DWKDX(IZ))/2/R + +D(NE+1,NE,IZ,2,6)*SF(I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ))/R + +D(NE+1,NE,IZ,2,2)*GSF(2,J)*SF(I)*(DWDY(IZ)+2*DWKDY(IZ))/2/R + +DB(NE,IZ,NE+1,1,3)*SF(I)*GSF(1,J)*(DWDX(IZ)+2*DWKDX(IZ))/2 + +DB(NE,IZ,NE+1,2,3)*GSF(2,J)*SF(I)*(DWDY(IZ)+2*DWKDY(IZ))/2 + +DB(NE,IZ,NE+1,3,6)*SF(I)*GSF(1,J)*(DWDY(IZ)+DWKDY(IZ)) + + +D(NE+1,NE,IZ,1,2)*GSF(1,I)*SF(J)*(DWDX(IZ)+DWKDX(IZ))/R + +D(NE+1,NE,IZ,2,2)*GSF(2,I)*SF(J)*(DWDY(IZ)+DWKDY(IZ))/R + +DB(NE+1,NE,IZ,1,3)*GSF(1,I)*SF(J)*(DWDX(IZ)+DWKDX(IZ)) + +DB(NE+1,NE,IZ,2,3)*GSF(2,I)*SF(J)*(DWDY(IZ)+DWKDY(IZ)) + +DB(NE+1,NE,IZ,3,6)*SF(J)*(GSF(1,I)*(DWDY(IZ)+DWKDY(IZ)) + +GSF(2,I)*(DWDX(IZ)+DWKDX(IZ))) + +D(NE+1,NE,IZ,2,6)*SF(J)*(GSF(1,I)*(DWDY(IZ)+DWKDY(IZ)) + + GSF(2,I)*(DWDX(IZ)+DWKDX(IZ)))/R IZ = IZ+1 200 CONTINUE END IF C C.... DOUBLE NON-LINEAR TERMS OF K33 FOLLOW C IF(NE.EQ.1) THEN NC = 2 JE = NE KE = NE ELSE NC = 3 JE = NE-1 KE = NE-1 END IF DO 220 IC = 1,NC NU = 2 IF(NE.NE.1.AND.IC.EQ.2) NU = 3 DO 225 ID = 1,NU IF(JE.GT.10.OR.KE.GT.10) THEN JE = NE KE = NE END IF TERM1 = 0.5*(DWDX(JE)+DWKDX(JE))*GSF(1,I)* + (F(NE,NE,JE,KE,1,1)*GSF(1,J)*(DWDX(KE)+2*DWKDX(KE))+ + F(NE,NE,JE,KE,1,2)*GSF(2,J)*(DWDY(KE)+2*DWKDY(KE))+ + F(NE,NE,JE,KE,1,6)*(GSF(1,J)*(DWDY(KE)+2*DWKDY(KE))+ + GSF(2,J)*(DWDX(KE)+2*DWKDX(KE)))) TERM2 = 0.5*(DWDY(JE)+DWDY(JE))*GSF(2,I)* + (F(NE,NE,JE,KE,1,2)*GSF(1,J)*(DWDX(KE)+2*DWKDX(KE))+ + F(NE,NE,JE,KE,2,2)*GSF(2,J)*(DWDY(KE)+2*DWKDY(KE))+ + F(NE,NE,JE,KE,2,6)*(GSF(1,J)*(DWDY(KE)+2*DWKDY(KE))+ + GSF(2,J)*(DWDX(KE)+2*DWKDX(KE)))) EXT = 0.5*(GSF(1,I)*(DWDY(JE)+DWKDY(JE))+ + GSF(2,I)*(DWDX(JE)+DWKDX(JE))) TERM3 = EXT*(F(NE,NE,JE,KE,1,6)*GSF(1,J)*(DWDX(KE)+2*DWKDX(KE))+ + F(NE,NE,JE,KE,2,6)*GSF(2,J)*(DWDY(KE)+2*DWKDY(KE))+ + F(NE,NE,JE,KE,6,6)*(GSF(1,J)*(DWDY(KE)+2*DWKDY(KE))+ + GSF(2,J)*(DWDX(KE)+2*DWKDX(KE)))) CON555 = CON555 + TERM1 + TERM2 + TERM3 225 KE = KE+1 KE = KE-1 JE = JE+1 220 CONTINUE IF(IN.LT.JN) THEN JE = NE KE = NE DO 230 IC = 1,2 DO 235 ID = 1,2 IF(JE.GT.10.OR.KE.GT.10) THEN JE = NE KE = NE END IF TERM4 = 0.5*(DWDX(JE)+DWKDX(JE))*GSF(1,I)* + (F(NE,NE+1,JE,KE,1,1)*GSF(1,J)*(DWDX(KE)+2*DWKDX(KE))+ + F(NE,NE+1,JE,KE,1,2)*GSF(2,J)*(DWDY(KE)+2*DWKDY(KE))+ + F(NE,NE+1,JE,KE,1,6)*(GSF(1,J)*(DWDY(KE)+2*DWKDY(KE))+ + GSF(2,J)*(DWDX(KE)+2*DWKDX(KE)))) TERM5 = 0.5*(DWDY(JE)+DWDY(JE))*GSF(2,I)* + (F(NE,NE+1,JE,KE,1,2)*GSF(1,J)*(DWDX(KE)+2*DWKDX(KE))+ + F(NE,NE+1,JE,KE,2,2)*GSF(2,J)*(DWDY(KE)+2*DWKDY(KE))+ + F(NE,NE+1,JE,KE,2,6)*(GSF(1,J)*(DWDY(KE)+2*DWKDY(KE))+ + GSF(2,J)*(DWDX(KE)+2*DWKDX(KE)))) EXT = 0.5*(GSF(1,I)*(DWDY(JE)+DWKDY(JE))+ + GSF(2,I)*(DWDX(JE)+DWKDX(JE))) TERM6 = EXT*(F(NE,NE+1,JE,KE,1,6)*GSF(1,J)*(DWDX(KE)+2*DWKDX(KE))+ + F(NE,NE+1,JE,KE,2,6)*GSF(2,J)*(DWDY(KE)+2*DWKDY(KE))+ + F(NE,NE+1,JE,KE,6,6)*(GSF(1,J)*(DWDY(KE)+2*DWKDY(KE)) + +GSF(2,J)*(DWDX(KE)+2*DWKDX(KE)))) CON666 = CON666 + TERM4 + TERM5 + TERM6 235 KE = KE+1 KE = NE JE = JE+1 230 CONTINUE END IF C C.... Contributions from 1/2[K33(Delta-GAMMA)] (Alpha,Gamma,Beta) C.... Follow C ME = NE-1 DO 192 IC=1,3 IF(ME.EQ.0) ME = NE CON55M = CON55M+ + D(NE,ME,NE,1,2)*GSF(1,J)*SF(I)*DWDX(ME)/2/R + +D(NE,ME,NE,2,6)*SF(I)*GSF(2,J)*DWDY(ME)/R + +D(NE,ME,NE,2,2)*GSF(2,J)*SF(I)*DWDY(ME)/2/R + +DB(ME,NE,NE,1,3)*SF(I)*GSF(1,J)*DWDX(ME)/2 + +DB(ME,NE,NE,2,3)*GSF(2,J)*SF(I)*DWDY(ME)/2 + +DB(ME,NE,NE,3,6)*SF(I)*GSF(1,J)*DWDY(ME) + +D(NE,ME,NE,1,2)*GSF(1,I)*SF(J)*DWDX(ME)/R + +D(NE,ME,NE,2,2)*GSF(2,I)*SF(J)*DWDY(ME)/R + +DB(NE,ME,NE,1,3)*GSF(1,I)*SF(J)*DWDX(ME) + +DB(NE,ME,NE,2,3)*GSF(2,I)*SF(J)*DWDY(ME) + +DB(NE,ME,NE,3,6)*SF(J)*(GSF(1,I)*DWDY(ME)+GSF(2,I)*DWDX(ME)) + +D(NE,ME,NE,2,6)*SF(J) *(GSF(1,I)*DWDY(ME)+GSF(2,I)*DWDX(ME))/R ME = ME+1 192 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 202 IC=1,2 CON66M = CON66M + D(NE,IZ,NE+1,1,2)*GSF(1,J)*SF(I)*DWDX(IZ)/2/R + + D(NE,IZ,NE+1,2,6)*SF(I)*GSF(2,J)*DWDY(IZ)/R + + D(NE,IZ,NE+1,2,2)*GSF(2,J)*SF(I)*DWDY(IZ)/2/R + + DB(IZ,NE+1,NE,1,3)*SF(I)*GSF(1,J)*DWDX(IZ)/2 + + DB(IZ,NE+1,NE,2,3)*GSF(2,J)*SF(I)*DWDY(IZ)/2 + + DB(IZ,NE+1,NE,3,6)*SF(I)*GSF(1,J)*DWDY(IZ) + + D(NE,IZ,NE+1,1,2)*GSF(1,I)*SF(J)*DWDX(IZ)/R + + D(NE,IZ,NE+1,2,2)*GSF(2,I)*SF(J)*DWDY(IZ)/R + + DB(NE,IZ,NE+1,1,3)*GSF(1,I)*SF(J)*DWDX(IZ) + + DB(NE,IZ,NE+1,2,3)*GSF(2,I)*SF(J)*DWDY(IZ) + +DB(NE,IZ,NE+1,3,6)*SF(J)*(GSF(1,I)*DWDY(IZ)+GSF(2,I)*DWDX(IZ)) + +D(NE,IZ,NE+1,2,6)*SF(J) *(GSF(1,I)*DWDY(IZ)+GSF(2,I)*DWDX(IZ))/R CON66J = CON66J + D(IZ,NE,NE+1,1,2)*GSF(1,J)*SF(I)*DWDX(IZ)/2/R + + D(IZ,NE,NE+1,2,6)*SF(I)*GSF(2,J)*DWDY(IZ)/R + + D(IZ,NE,NE+1,2,2)*GSF(2,J)*SF(I)*DWDY(IZ)/2/R + + DB(NE,NE+1,IZ,1,3)*SF(I)*GSF(1,J)*DWDX(IZ)/2 + + DB(NE,NE+1,IZ,2,3)*GSF(2,J)*SF(I)*DWDY(IZ)/2 + + DB(NE,NE+1,IZ,3,6)*SF(I)*GSF(1,J)*DWDY(IZ) + + D(IZ,NE,NE+1,1,2)*GSF(1,I)*SF(J)*DWDX(IZ)/R + + D(IZ,NE,NE+1,2,2)*GSF(2,I)*SF(J)*DWDY(IZ)/R + + DB(IZ,NE,IZ+1,1,3)*GSF(1,I)*SF(J)*DWDX(IZ) + + DB(IZ,NE,IZ+1,2,3)*GSF(2,I)*SF(J)*DWDY(IZ) + +DB(IZ,NE,NE+1,3,6)*SF(J)*(GSF(1,I)*DWDY(IZ)+GSF(2,I)*DWDX(IZ)) + +D(IZ,NE,NE+1,2,6)*SF(J)*(GSF(1,I)*DWDY(IZ)+GSF(2,I)*DWDX(IZ))/R IZ = IZ+1 202 CONTINUE END IF C C.... Double nonlinear contributions from 1/2[K33(Dq + Dbarq)] C follow IF(NE.EQ.1) THEN NC = 2 JE = NE KE = NE ELSE NC = 3 JE = NE-1 KE = NE-1 END IF DO 222 IC = 1,NC NU = 2 IF(NE.NE.1.AND.IC.EQ.2) NU = 3 DO 227 ID = 1,NU IF(JE.GT.10.OR.KE.GT.10) THEN JE = NE KE = NE END IF TERM1 = 0.5*(DWDX(KE)+DWKDX(KE))*GSF(1,I)* + (F(NE,JE,NE,KE,1,1)*GSF(1,J)*DWDX(JE)+ + F(NE,JE,NE,KE,1,2)*GSF(2,J)*DWDY(JE)+ + F(NE,JE,NE,KE,1,6)*(GSF(1,J)*DWDY(JE)+GSF(2,J)*DWDX(JE))) TERM2 = 0.5*(DWDY(KE)+DWKDY(KE))*GSF(2,I)* + (F(NE,JE,NE,KE,1,2)*GSF(1,J)*DWDX(JE)+ + F(NE,JE,NE,KE,2,2)*GSF(2,J)*DWDY(JE)+ + F(NE,JE,NE,KE,2,6)*(GSF(1,J)*DWDY(JE)+GSF(2,J)*DWDX(JE))) EXT = 0.5*(GSF(1,I)*DWDY(JE)+GSF(2,I)*DWDX(JE)) TERM3 = EXT*(F(NE,JE,NE,KE,1,6)*GSF(1,J)*(DWDX(KE)+DWKDX(KE))+ + F(NE,JE,NE,KE,2,6)*GSF(2,J)*(DWDY(KE)+DWKDY(KE))+ + F(NE,JE,NE,KE,6,6)*(GSF(1,J)*(DWDY(KE)+DWKDY(KE))+ + GSF(2,J)*(DWDX(KE)+DWKDX(KE)))) CON55T = CON55T + TERM1 + TERM2 + TERM3 227 KE = KE+1 KE = KE-1 JE = JE+1 222 CONTINUE IF(IN.LT.JN) THEN JE = NE KE = NE DO 232 IC = 1,2 DO 237 ID = 1,2 IF(JE.GT.10.OR.KE.GT.10) THEN JE = NE KE = NE END IF TERM4 = 0.5*(DWDX(KE)+DWKDX(KE))*GSF(1,I)* + (F(NE,JE,NE+1,KE,1,1)*GSF(1,J)*DWDX(JE)+ + F(NE,JE,NE+1,KE,1,2)*GSF(2,J)*DWDY(JE) + +F(NE,JE,NE+1,KE,1,6)*(GSF(1,J)*DWDY(JE)+ + GSF(2,J)*DWDX(JE))) TERM5 = 0.5*(DWDY(JE)+DWDY(JE))*GSF(2,I)* + (F(NE,JE,NE+1,KE,1,2)*GSF(1,J)*DWDX(JE)+ + F(NE,JE,NE+1,KE,2,2)*GSF(2,J)*DWDY(JE)+ + F(NE,JE,NE+1,KE,2,6)*(GSF(1,J)*DWDY(JE)+ + GSF(2,J)*DWDX(KE))) EXT = 0.5*(GSF(1,I)*DWDY(JE)+GSF(2,I)*DWDX(JE)) TERM6 = EXT*(F(NE,JE,NE+1,KE,1,6)*GSF(1,J)*(DWDX(KE)+DWKDX(KE))+ + F(NE,JE,NE+1,KE,2,6)*GSF(2,J)*(DWDY(KE)+DWKDY(KE))+ + F(NE,JE,NE+1,KE,6,6)*(GSF(1,J)*(DWDY(KE)+DWKDY(KE)) + + GSF(2,J)*(DWDX(KE)+DWKDX(KE)))) CON66T = CON66T + TERM4 + TERM5 + TERM6 237 KE = KE+1 KE = NE JE = JE+1 232 CONTINUE END IF C C....MORE DOUBLE NONLINEAR CONTRIBUTIONS FROM 1/2[K33(DRHO/2+ DBARRHO)] C FOLOW. ( Alpha, Gamma, Rho, Beta) IF(NE.EQ.1) THEN NC = 2 JE = NE KE = NE ELSE NC = 3 JE = NE-1 KE = NE-1 END IF DO 224 IC = 1,NC NU = 2 IF(NE.NE.1.AND.IC.EQ.2) NU = 3 DO 229 ID = 1,NU IF(JE.GT.10.OR.KE.GT.10) THEN JE = NE KE = NE END IF TERM1 = 0.5*(DWDX(KE)+2*DWKDX(KE))*GSF(1,I)* + (F(NE,JE,KE,NE,1,1)*GSF(1,J)*DWDX(JE)+ + F(NE,JE,KE,NE,1,2)*GSF(2,J)*DWDY(JE)+ + F(NE,JE,KE,NE,1,6)*(GSF(1,J)*DWDY(JE)+GSF(2,J)*DWDX(JE))) TERM2 = 0.5*(DWDY(KE)+2*DWKDY(KE))*GSF(2,I)* + (F(NE,JE,KE,NE,1,2)*GSF(1,J)*DWDX(JE)+ + F(NE,JE,KE,NE,2,2)*GSF(2,J)*DWDY(JE)+ + F(NE,JE,KE,NE,2,6)*(GSF(1,J)*DWDY(JE)+GSF(2,J)*DWDX(JE))) EXT = 0.5*(GSF(1,I)*DWDY(JE)+GSF(2,I)*DWDX(JE)) TERM3 = EXT*(F(NE,JE,KE,NE,1,6)*GSF(1,J)*(DWDX(KE)+2*DWKDX(KE))+ + F(NE,JE,KE,NE,2,6)*GSF(2,J)*(DWDY(KE)+2*DWKDY(KE))+ + F(NE,JE,KE,NE,6,6)*(GSF(1,J)*(DWDY(KE)+2*DWKDY(KE))+ + GSF(2,J)*(DWDX(KE)+2*DWKDX(KE)))) CON55S = CON55S + TERM1 + TERM2 + TERM3 229 KE = KE+1 KE = KE-1 JE = JE+1 224 CONTINUE IF(IN.LT.JN) THEN JE = NE KE = NE DO 234 IC = 1,2 DO 239 ID = 1,2 IF(JE.GT.10.OR.KE.GT.10) THEN JE = NE KE = NE END IF TERM4 = 0.5*(DWDX(KE)+2*DWKDX(KE))*GSF(1,I)* + (F(NE,JE,KE,NE+1,1,1)*GSF(1,J)*DWDX(JE)+ + F(NE,JE,KE,NE+1,1,2)*GSF(2,J)*DWDY(JE)+ + F(NE,JE,KE,NE+1,1,6)*(GSF(1,J)*DWDY(JE)+ + GSF(2,J)*DWDX(JE))) TERM5 = 0.5*(DWDY(JE)+2*DWDY(JE))*GSF(2,I)* + (F(NE,JE,KE,NE+1,1,2)*GSF(1,J)*DWDX(JE)+ + F(NE,JE,KE,NE+1,2,2)*GSF(2,J)*DWDY(JE)+ + F(NE,JE,KE,NE+1,2,6)*(GSF(1,J)*DWDY(JE)+ + GSF(2,J)*DWDX(KE))) EXT = 0.5*(GSF(1,I)*DWDY(JE)+GSF(2,I)*DWDX(JE)) TERM6 = EXT*(F(NE,JE,KE,NE+1,1,6)*GSF(1,J)*(DWDX(KE)+2*DWKDX(KE))+ + F(NE,JE,KE,NE+1,2,6)*GSF(2,J)*(DWDY(KE)+2*DWKDY(KE))+ + F(NE,JE,KE,NE+1,6,6)*(GSF(1,J)*(DWDY(KE)+2*DWKDY(KE))+ + GSF(2,J)*(DWDX(KE)+2*DWKDX(KE)))) CON66S = CON66S + TERM4 + TERM5 + TERM6 239 KE = KE+1 KE = NE JE = JE+1 234 CONTINUE END IF C C.... NON-LINEAR CONTRIBUTIONS FROM (DK31/DW)*UJ FOLLOW C ME = NE-1 DO 552 IC = 1,3 IF (ME.EQ.0) ME = NE CON311 = CON311 + + D(NE,NE,ME,1,1)*GSF(1,J)*GSF(1,I)*DUDX(ME) + +D(NE,NE,ME,1,2)*GSF(1,J)*GSF(2,I)*DUDY(ME) + +D(NE,NE,ME,1,6)*GSF(2,J)*GSF(1,I)*DUDX(ME) + +D(NE,NE,ME,2,6)*GSF(2,J)*GSF(2,I)*DUDY(ME) + +D(NE,NE,ME,1,6)*GSF(1,J)*(GSF(2,I)*DUDX(ME)+GSF(1,I)*DUDY(ME)) + +D(NE,NE,ME,6,6)*GSF(2,J)*(GSF(2,I)*DUDX(ME)+GSF(1,I)*DUDY(ME)) ME = ME+1 552 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 562 IC = 1,2 CON312 = CON312 + + D(NE,NE+1,IZ,1,1)*GSF(1,J)*GSF(1,I)*DUDX(IZ) + + D(NE,NE+1,IZ,1,2)*GSF(1,J)*GSF(2,I)*DUDY(IZ) + + D(NE,NE+1,IZ,1,6)*GSF(2,J)*GSF(1,I)*DUDX(IZ) + + D(NE,NE+1,IZ,2,6)*GSF(2,J)*GSF(2,I)*DUDY(IZ) + + D(NE,NE+1,IZ,1,6)*GSF(1,J)*(GSF(2,I)*DUDX(IZ)+GSF(1,I)*DUDX(IZ))+ + D(NE,NE+1,IZ,6,6)*GSF(2,J)*(GSF(2,I)*DUDX(IZ)+GSF(1,I)*DUDY(IZ)) IZ = IZ+1 562 CONTINUE END IF C C.... NON-LINEAR TERMS CONTRIBUTIONS FROM (DK32/DW)*VJ FOLLOW C ME = NE-1 DO 570 IC=1,3 IF(ME.EQ.0) ME = NE CON321 = CON321 + + D(NE,NE,ME,1,2)*GSF(2,J)*GSF(1,I)*DVDX(ME) + +D(NE,NE,ME,2,2)*GSF(2,J)*GSF(2,I)*DVDY(ME) + +D(NE,NE,ME,1,6)*GSF(1,J)*GSF(1,I)*DVDX(ME) + +D(NE,NE,ME,2,6)*GSF(1,J)*GSF(2,I)*DVDY(ME) + +D(NE,NE,ME,2,6)*GSF(2,J)*(GSF(1,I)*DVDY(ME)+GSF(2,I)*DVDX(ME)) + + D(NE,NE,ME,6,6)*GSF(2,J)*(GSF(1,I)*DVDY(ME)+GSF(2,I)*DVDX(ME)) 574 ME = ME+1 570 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 580 IC=1,2 CON322 = CON322 + + D(NE,NE+1,IZ,1,2)*GSF(2,J)*GSF(1,I)*DVDX(IZ) + +D(NE,NE+1,IZ,2,2)*GSF(2,J)*GSF(2,I)*DVDY(IZ) + +D(NE,NE+1,IZ,1,6)*GSF(1,J)*GSF(1,I)*DVDX(IZ) + +D(NE,NE+1,IZ,2,6)*GSF(1,J)*GSF(2,I)*DVDY(IZ)+ + D(NE,NE+1,IZ,2,6)*GSF(2,J)*(GSF(1,I)*DVDY(IZ)+GSF(2,I)*DVDX(IZ))+ + D(NE,NE+1,IZ,6,6)*GSF(2,J)*(GSF(1,I)*DVDY(IZ)+GSF(2,I)*DVDX(IZ)) IZ = IZ+1 580 CONTINUE END IF ST(IQ,IP) = ST(IQ,IP) +(CON5+CON55+CON555)*CNST TANK(IQ,IP) = TANK(IQ,IP)+(CON5+CON311+CON321+CON55+ + CON555+CON55M+CON55T+CON55S)*CNST 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) +(CON6+CON66I+CON666)*CNST TANK(IQ,IP+1)= TANK(IQ,IP+1)+(CON6+CON312+CON322+CON66+ + CON666+CON66M+CON55T+CON66S)*CNST TANK(IQ+1,IP)= TANK(IQ+1,IP)+(CON6+CON312+CON322+CON66I+ + CON666+CON66M+CON55T+CON66S)*CNST END IF C_______________________________________________________________________ C C This segment of the program models imperfections in the C cylindrical shell. It is activated for IMPFCT = 1 only. C C_______________________________________________________________________ IF(IMPFCT.EQ.1) THEN C_______________________________________________________________________ C C.....Calculate [K14] (For Imperfection Sensitivity Analysis only) C_______________________________________________________________________ IQ = I1 IP = IN + 3*NIF ME = NE-1 DO 650 IC=1,3 IF(ME.EQ.0) ME = NE CON14 = CON14 + + D(NE,NE,ME,1,1)*GSF(1,I)*GSF(1,J)*DWDX(NE) + +D(NE,NE,ME,1,6)*GSF(1,I)*GSF(2,J)*DWDX(NE) + +D(NE,NE,ME,1,2)*GSF(1,I)*GSF(2,J)*DWDY(NE) + +D(NE,NE,ME,1,6)*GSF(2,I)*GSF(1,J)*DWDX(NE) + +D(NE,NE,ME,2,6)*GSF(2,I)*GSF(2,J)*DWDY(NE) + +D(NE,NE,ME,6,6)*GSF(2,I)*GSF(2,J)*DWDX(NE)/2 ME = ME+1 650 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 660 IC=1,2 CON141 = CON141 + + D(NE,NE+1,IZ,1,1)*GSF(1,I)*GSF(1,J)*DWDX(NE+1) + +D(NE,NE+1,IZ,1,6)*GSF(1,I)*GSF(2,J)*DWDX(NE+1) + +D(NE,NE+1,IZ,1,2)*GSF(1,I)*GSF(2,J)*DWDY(NE+1) + +D(NE,NE+1,IZ,1,6)*GSF(2,I)*GSF(1,J)*DWDX(NE+1) + +D(NE,NE+1,IZ,2,6)*GSF(2,I)*GSF(2,J)*DWDY(NE+1) + +D(NE,NE+1,IZ,6,6)*GSF(2,I)*GSF(2,J)*DWDX(NE+1)/2 CON14I= CON14I + + D(NE+1,NE,IZ,1,1)*GSF(1,I)*GSF(1,J)*DWDX(NE) + +D(NE+1,NE,IZ,1,6)*GSF(1,I)*GSF(2,J)*DWDX(NE) + +D(NE+1,NE,IZ,1,2)*GSF(1,I)*GSF(2,J)*DWDY(NE) + +D(NE+1,NE,IZ,1,6)*GSF(2,I)*GSF(1,J)*DWDX(NE) + +D(NE+1,NE,IZ,2,6)*GSF(2,I)*GSF(2,J)*DWDY(NE) + +D(NE+1,NE,IZ,6,6)*GSF(2,I)*GSF(2,J)*DWDX(NE)/2 IZ = IZ+1 660 CONTINUE END IF ST(IQ,IP) = ST(IQ,IP) + CON14*CNST TANK(IQ,IP) = TANK(IQ,IP)+ 2*CON14*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1)+CON141*CNST ST(IQ+1,IP) = ST(IQ+1,IP)+ CON14I*CNST TANK(IQ,IP+1)= TANK(IQ,IP+1)+2*CON141*CNST TANK(IQ+1,IP)= TANK(IQ+1,IP)+ 2*CON14I*CNST END IF C_______________________________________________________________________ C C.....Calculate [K24] (For Imperfection Sensitivity Analysis only) C______________________________________________________________________ IQ = I1 + NIF IP = IN + 3*NIF ME = NE-1 DO 670 IC=1,3 IF(ME.EQ.0) ME = NE CON24 = CON24+ + D(NE,NE,ME,1,2)*GSF(2,I)*GSF(1,J)*DWDX(NE) + +D(NE,NE,ME,2,6)*GSF(2,I)*GSF(2,J)*DWDX(NE) + +D(NE,NE,ME,2,2)*GSF(2,I)*GSF(2,J)*DWDY(NE) + +D(NE,NE,ME,1,6)*GSF(1,I)*GSF(1,J)*DWDX(NE) + +D(NE,NE,ME,2,6)*GSF(1,I)*GSF(2,J)*DWDY(NE) + +D(NE,NE,ME,6,6)*GSF(1,I)*GSF(2,J)*DWDX(NE)/2 ME = ME+1 670 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 680 IC=1,2 CON242 = CON242 + + D(NE,NE+1,IZ,1,2)*GSF(2,I)*GSF(1,J)*DWDX(NE+1) + +D(NE,NE+1,IZ,2,6)*GSF(2,I)*GSF(2,J)*DWDX(NE+1) + +D(NE,NE+1,IZ,2,2)*GSF(2,I)*GSF(2,J)*DWDY(NE+1) + +D(NE,NE+1,IZ,1,6)*GSF(1,I)*GSF(1,J)*DWDX(NE+1) + +D(NE,NE+1,IZ,2,6)*GSF(1,I)*GSF(2,J)*DWDY(NE+1) + +D(NE,NE+1,IZ,6,6)*GSF(1,I)*GSF(2,J)*DWDX(NE+1)/2 IZ = IZ+1 680 CONTINUE END IF ST(IQ,IP) = ST(IQ,IP) + CON24*CNST TANK(IQ,IP) = TANK(IQ,IP)+ 2*CON24*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1) + CON242*CNST ST(IQ+1,IP) = ST(IQ+1,IP) + CON242*CNST TANK(IQ,IP+1) = TANK(IQ,IP+1)+ 2*CON242*CNST TANK(IQ+1,IP) = TANK(IQ+1,IP)+ 2*CON242*CNST END IF C_______________________________________________________________________ C C.....Calculate [K34] (For Imperfection Sensitivity Analysis only) C______________________________________________________________________ IQ = I1 + 2*NIF IP = IN + 3*NIF ME = NE-1 DO 690 IC=1,3 IF(ME.EQ.0) ME = NE CON34 = CON34 + + D(NE,NE,ME,1,2)*GSF(1,J)*SF(I)*DWDX(NE)/R + +D(NE,NE,ME,2,6)*SF(I)*GSF(2,J)*DWDX(NE)/R + +D(NE,NE,ME,2,2)*GSF(2,J)*SF(I)*DWDY(NE)/R + +DB(NE,ME,NE,1,3)*SF(I)*GSF(1,J)*DWDX(NE) + +DB(NE,ME,NE,2,3)*SF(I)*GSF(2,J)*DWDY(NE) + +DB(NE,ME,NE,3,6)*SF(I)*GSF(2,J)*DWDX(NE) CON34A= CON34A+ + D(NE,ME,NE,1,2)*GSF(1,J)*SF(I)*DWDX(NE)/R + +D(NE,ME,NE,2,6)*SF(I)*GSF(2,J)*DWDX(NE)/R + +D(NE,ME,NE,2,2)*GSF(2,J)*SF(I)*DWDY(NE)/R + +DB(NE,NE,ME,1,3)*SF(I)*GSF(1,J)*DWDX(NE) + +DB(NE,NE,ME,2,3)*SF(I)*GSF(2,J)*DWDY(NE) + +DB(NE,NE,ME,3,6)*SF(I)*GSF(2,J)*DWDX(NE) ME = ME+1 690 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 700 IC=1,2 CON341= CON341 + + D(NE,NE+1,IZ,1,2)*GSF(1,J)*SF(I)*DWDX(NE+1)/R + +D(NE,NE+1,IZ,2,6)*SF(I)*GSF(2,J)*DWDX(NE+1)/R + +D(NE,NE+1,IZ,2,2)*GSF(2,J)*SF(I)*DWDY(NE+1)/R + +DB(NE+1,IZ,NE,1,3)*SF(I)*GSF(1,J)*DWDX(NE+1) + +DB(NE+1,IZ,NE,2,3)*SF(I)*GSF(2,J)*DWDY(NE+1) + +DB(NE+1,IZ,NE,3,6)*SF(I)*GSF(2,J)*DWDX(NE+1) CON34B= CON34B + + D(NE,NE+1,IZ,1,2)*GSF(1,J)*SF(I)*DWDX(NE+1)/R + +D(NE,NE+1,IZ,2,6)*SF(I)*GSF(2,J)*DWDX(NE+1)/R + +D(NE,NE+1,IZ,2,2)*GSF(2,J)*SF(I)*DWDY(NE+1)/R + +DB(NE+1,NE,IZ,1,3)*SF(I)*GSF(1,J)*DWDX(NE+1) + +DB(NE+1,NE,IZ,2,3)*SF(I)*GSF(2,J)*DWDY(NE+1) + +DB(NE+1,NE,IZ,3,6)*SF(I)*GSF(2,J)*DWDX(NE+1) CON34I = CON34I + + D(NE+1,NE,IZ,1,2)*GSF(1,J)*SF(I)*DWDX(NE)/R + +D(NE+1,NE,IZ,2,6)*SF(I)*GSF(2,J)*DWDX(NE)/R + +D(NE+1,NE,IZ,2,2)*GSF(2,J)*SF(I)*DWDY(NE)/R + +DB(NE,IZ,NE+1,1,3)*SF(I)*GSF(1,J)*DWDX(NE) + +DB(NE,IZ,NE+1,2,3)*SF(I)*GSF(2,J)*DWDY(NE) + +DB(NE,IZ,NE+1,3,6)*SF(I)*GSF(2,J)*DWDX(NE) CON34C = CON34C + + D(NE+1,NE,IZ,1,2)*GSF(1,J)*SF(I)*DWDX(NE)/R + +D(NE+1,NE,IZ,2,6)*SF(I)*GSF(2,J)*DWDX(NE)/R + +D(NE+1,NE,IZ,2,2)*GSF(2,J)*SF(I)*DWDY(NE)/R + +DB(NE,NE+1,IZ,1,3)*SF(I)*GSF(1,J)*DWDX(NE) + +DB(NE,NE+1,IZ,2,3)*SF(I)*GSF(2,J)*DWDY(NE) + +DB(NE,NE+1,IZ,3,6)*SF(I)*GSF(2,J)*DWDX(NE) IZ = IZ+1 700 CONTINUE END IF C C.... Double non-linear terms of [K34] ...... C IF(NE.EQ.1) THEN NC = 2 JE = NE KE = NE ELSE NC = 3 JE = NE-1 KE = NE-1 END IF DO 720 IC = 1,NC NU = 2 IF(NE.NE.1.AND.IC.EQ.2) NU = 3 DO 725 ID = 1,NU TERM1 = GSF(1,I)*DWKDX(NE)* + (F(NE,NE,JE,KE,1,1)*GSF(1,J)*DWDX(JE)+ + F(NE,NE,JE,KE,1,2)*GSF(2,J)*DWDY(JE)+ + F(NE,NE,JE,KE,1,6)*(GSF(1,J)*DWDY(JE)+GSF(2,J)*DWDX(JE))) TERM2 = GSF(2,I)*DWKDY(NE)* + (F(NE,NE,JE,KE,1,2)*GSF(1,J)*DWDX(JE)+ + F(NE,NE,JE,KE,2,2)*GSF(2,J)*DWDY(JE)+ + F(NE,NE,JE,KE,2,6)*(GSF(1,J)*DWDY(JE)+GSF(2,J)*DWDX(JE))) EXT = (GSF(1,I)*DWKDY(NE)+GSF(2,I)*DWKDX(NE)) TERM3 = EXT*(F(NE,NE,JE,KE,1,6)*GSF(1,J)*DWDX(JE)+ + F(NE,NE,JE,KE,2,6)*GSF(2,J)*DWDY(JE)+ + F(NE,NE,JE,KE,6,6)*(GSF(1,J)*DWDY(JE)+ + GSF(2,J)*DWDX(JE))) CON344 = CON344 + TERM1 + TERM2 +TERM3 725 KE = KE+1 KE = KE-1 JE = JE+1 720 CONTINUE IF(IN.LT.JN) THEN JE = NE KE = NE DO 730 IC = 1,2 DO 735 ID = 1,2 TERM4 = GSF(1,I)*DWKDX(NE+1)* + (F(NE,NE+1,JE,KE,1,1)*GSF(1,J)*DWDX(JE)+ + F(NE,NE+1,JE,KE,1,2)*GSF(2,J)*DWDY(JE)+ + F(NE,NE+1,JE,KE,1,6)*(GSF(1,J)*DWDY(JE)+GSF(2,J)*DWDX(JE))) TERM5 = GSF(2,I)*DWKDY(NE)* + (F(NE,NE+1,JE,KE,1,2)*GSF(1,J)*DWDX(JE)+ + F(NE,NE+1,JE,KE,2,2)*GSF(2,J)*DWDY(JE)+ + F(NE,NE+1,JE,KE,2,6)*(GSF(1,J)*DWDY(JE)+GSF(2,J)*DWDX(JE))) EXT = (GSF(1,I)*DWKDY(NE+1)+ GSF(2,I)*DWKDX(NE+1)) TERM6 = EXT*(F(NE,NE+1,JE,KE,1,6)*GSF(1,J)*DWDX(JE)+ + F(NE,NE+1,JE,KE,2,6)*GSF(2,J)*DWDY(JE)+ + F(NE,NE+1,JE,KE,6,6)*(GSF(1,J)*DWDY(JE)+ + GSF(2,J)*DWDX(JE))) CON355 = CON355+TERM4 + TERM5 +TERM6 735 KE = KE+1 KE = NE JE = JE+1 730 CONTINUE END IF ST(IQ,IP) = ST(IQ,IP) + (CON34+CON344)*CNST TANK(IQ,IP) = TANK(IQ,IP)+ (CON34+CON344+CON34A+CON344)*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1) + (CON341+CON355)*CNST ST(IQ+1,IP) = ST(IQ+1,IP) + (CON34I+CON355)*CNST TANK(IQ,IP+1) = TANK(IQ,IP+1)+ (CON341+CON355+CON34B+CON355)*CNST TANK(IQ+1,IP) = TANK(IQ+1,IP)+ (CON34I+CON355+CON34C+CON355)*CNST END IF END IF C_______________________________________________________________________ C C Imperfection Sensitivity Analysis Ends Here.... C_______________________________________________________________________ I1 = I1+1 55 CONTINUE 60 JJ= NDF*J+1 70 II= NDF*I+1 100 CONTINUE C---------------------------------------------------------------------- C CARRY OUT REDUCED INTEGRATION FOR TRANSVERSE SHEAR TERMS C (Terms containing Qi4's and Qi5's, i = 4,5) C C---------------------------------------------------------------------- DO 600 NI =1,LGP DO 600 NJ =1,LGP XI = GAUSS(NI,LGP) ETA = GAUSS(NJ,LGP) CALL SHAPE2(NPE,XI,ETA,SF,GSF,DET,ELXY) CNST= DET*WT(NI,LGP)*WT(NJ,LGP) II = 1 DO 770 I = 1,NPE JJ = 1 DO 860 J = 1,NPE I1 = II JN=JJ+NPLY NE = 0 DO 655 IN = JJ,JN NE = NE+1 C------------------------------------------------------------------- C.....CALCULATE [K11] -- Shear terms C-------------------------------------------------------------------- ST(I1,IN) = ST(I1,IN) + DBAR(NE,NE,5,5)*SF(I)*SF(J)*CNST TANK(I1,IN) = ST(I1,IN) IF(IN.LT.JN) THEN ST(I1,IN+1) = ST(I1,IN+1)+ DBAR(NE,NE+1,5,5)*SF(I)*SF(J)*CNST ST(I1+1,IN) = ST(I1,IN+1) TANK(I1,IN+1) = ST(I1,IN+1) TANK(I1+1,IN) = ST(I1,IN+1) END IF C---------------------------------------------------------------------- C.....CALCULATE [K12] -- Shear terms C---------------------------------------------------------------------- ST(I1,IN+NIF) = ST(I1,IN+NIF)+(DBAR(NE,NE,4,5)*SF(I)*SF(J) + - ABAR(NE,NE,4,5)*SF(I)*SF(J)/R)*CNST TANK(I1,IN+NIF) = ST(I1,IN+NIF) IF(IN.LT.JN) THEN ST(I1,IN+NIF+1) = ST(I1,IN+NIF+1)+(DBAR(NE,NE+1,4,5)*SF(I)*SF(J) + - ABAR(NE+1,NE,4,5)*SF(I)*SF(J)/R)*CNST ST(I1+1,IN+NIF) = ST(I1,IN+NIF+1) TANK(I1+1,IN+NIF) = ST(I1,IN+NIF+1) TANK(I1,IN+NIF+1) = ST(I1,IN+NIF+1) END IF C---------------------------------------------------------------------- C.....CALCULATE [K21] - Shear terms C---------------------------------------------------------------------- ST(I1+NIF,IN) = ST(I1+NIF,IN)+(DBAR(NE,NE,4,5)*SF(I)*SF(J) + -ABAR(NE,NE,4,5)*SF(I)*SF(J)/R)*CNST TANK(I1+NIF,IN) = ST(I1+NIF,IN) IF(IN.LT.JN) THEN ST(I1+NIF,IN+1) = ST(I1+NIF,IN+1)+(DBAR(NE,NE+1,4,5)*SF(I)*SF(J) + -ABAR(NE+1,NE,4,5)*SF(I)*SF(J)/R)*CNST ST(I1+NIF+1,IN) = ST(I1+NIF,IN+1) TANK(I1+NIF,IN+1) = ST(I1+NIF,IN+1) TANK(I1+NIF+1,IN) = ST(I1+NIF,IN+1) END IF C------------------------------------------------------------------- C.....CALCULATE [K22] - Shear terms C-------------------------------------------------------------------- ST(I1+NIF,IN+NIF) = ST(I1+NIF,IN+NIF) + + (DBAR(NE,NE,4,4)*SF(I)*SF(J) + - ABAR(NE,NE,4,4)*SF(I)*SF(J)/R + - ABAR(NE,NE,4,4)*SF(I)*SF(J)/R + + A(NE,NE,4,4)*SF(I)*SF(J)/R/R)*CNST TANK(I1+NIF,IN+NIF) = ST(I1+NIF,IN+NIF) IF(IN.LT.JN) THEN ST(I1+NIF,IN+NIF+1) = ST(I1+NIF,IN+NIF+1) + +(DBAR(NE,NE+1,4,4)*SF(I)*SF(J) + - ABAR(NE+1,NE,4,4)*SF(I)*SF(J)/R + - ABAR(NE,NE+1,4,4)*SF(I)*SF(J)/R + + A(NE,NE+1,4,4)*SF(I)*SF(J)/R/R)*CNST ST(I1+NIF+1,IN+NIF) = ST(I1+NIF,IN+NIF+1) TANK(I1+NIF,IN+NIF+1) = ST(I1+NIF,IN+NIF+1) TANK(I1+NIF+1,IN+NIF) = ST(I1+NIF,IN+NIF+1) END IF C------------------------------------------------------------------- C CALCULATE [K13] - Shear terms C-------------------------------------------------------------------- IP = IN+2*NIF CON1 = ABAR(NE,NE,5,5)*SF(I)*GSF(1,J) + + ABAR(NE,NE,4,5)*SF(I)*GSF(2,J) IF(IN.LT.JN) THEN CON2 = ABAR(NE+1,NE,5,5)*SF(I)*GSF(1,J) + + ABAR(NE+1,NE,4,5)*SF(I)*GSF(2,J) CON2I= ABAR(NE,NE+1,5,5)*SF(I)*GSF(1,J) + +ABAR(NE,NE+1,4,5)*SF(I)*GSF(2,J) END IF ST(I1,IP) = ST(I1,IP)+CON1*CNST TANK(I1,IP) = TANK(I1,IP)+CON1*CNST IF(IN.LT.JN) THEN ST(I1,IP+1) = ST(I1,IP+1)+CON2*CNST ST(I1+1,IP) = ST(I1+1,IP)+CON2I*CNST TANK(I1,IP+1)= TANK(I1,IP+1)+CON2*CNST TANK(I1+1,IP)= TANK(I1+1,IP)+CON2I*CNST END IF C------------------------------------------------------------------- C CALCULATE [K31] -- Shear terms C-------------------------------------------------------------------- C C ....LINEAR TERMS FOLLOW C IP = IN IQ = I1+2*NIF CON1 = ABAR(NE,NE,5,5)*SF(J)*GSF(1,I) + + ABAR(NE,NE,4,5)*SF(J)*GSF(2,I) IF(IN.LT.JN) THEN CON2 = ABAR(NE,NE+1,5,5)*SF(J)*GSF(1,I) + + ABAR(NE,NE+1,4,5)*SF(J)*GSF(2,I) CON2I= ABAR(NE+1,NE,5,5)*SF(J)*GSF(1,I) + +ABAR(NE+1,NE,4,5)*SF(J)*GSF(2,I) END IF ST(IQ,IP) = ST(IQ,IP)+CON1*CNST TANK(IQ,IP) = TANK(IQ,IP)+CON1*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1)+CON2*CNST ST(IQ+1,IP) = ST(IQ+1,IP)+CON2I*CNST TANK(IQ,IP+1) = TANK(IQ,IP+1)+CON2*CNST TANK(IQ+1,IP) = TANK(IQ+1,IP)+CON2I*CNST END IF C------------------------------------------------------------------- C CALCULATE [K23] -- Shear terms C-------------------------------------------------------------------- IP = IN+2*NIF IQ = I1+NIF CON3 = ABAR(NE,NE,4,5)*SF(I)*GSF(1,J) + + ABAR(NE,NE,4,4)*SF(I)*GSF(2,J) + -A(NE,NE,4,5)*SF(I)*GSF(1,J)/R + -A(NE,NE,4,4)*SF(I)*GSF(2,J)/R IF(IN.LT.JN) THEN CON4 = ABAR(NE+1,NE,4,5)*SF(I)*GSF(1,J) + + ABAR(NE+1,NE,4,4)*SF(I)*GSF(2,J) + -A(NE,NE+1,4,5)*SF(I)*GSF(1,J)/R + -A(NE,NE+1,4,4)*SF(I)*GSF(2,J)/R CON4I= ABAR(NE,NE+1,4,5)*SF(I)*GSF(1,J) + + ABAR(NE,NE+1,4,4)*SF(I)*GSF(2,J) + - A(NE+1,NE,4,5)*SF(I)*GSF(1,J)/R + - A(NE+1,NE,4,4)*SF(I)*GSF(2,J)/R END IF ST(IQ,IP) = ST(IQ,IP) + CON3*CNST TANK(IQ,IP) = TANK(IQ,IP)+ CON3*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1)+ CON4*CNST ST(IQ+1,IP) = ST(IQ+1,IP)+ CON4I*CNST TANK(IQ,IP+1)= TANK(IQ,IP+1)+CON4*CNST TANK(IQ+1,IP)= TANK(IQ+1,IP)+CON4I*CNST END IF C------------------------------------------------------------------- C CALCULATE [K32] -- Shear terms C-------------------------------------------------------------------- IQ = I1+2*NIF IP = IN+NIF CON3 = ABAR(NE,NE,4,5)*SF(J)*GSF(1,I) + + ABAR(NE,NE,4,4)*SF(J)*GSF(2,I) + -A(NE,NE,4,5)*SF(J)*GSF(1,I)/R + -A(NE,NE,4,4)*SF(J)*GSF(2,I)/R IF(IN.LT.JN) THEN CON4 =ABAR(NE,NE+1,4,5)*SF(J)*GSF(1,I) + + ABAR(NE,NE+1,4,4)*SF(J)*GSF(2,I) + -A(NE,NE+1,4,5)*SF(J)*GSF(1,I)/R + -A(NE,NE+1,4,4)*SF(J)*GSF(2,I)/R CON4I=ABAR(NE+1,NE,4,5)*SF(J)*GSF(1,I) + + ABAR(NE+1,NE,4,4)*SF(J)*GSF(2,I) + -A(NE+1,NE,4,5)*SF(J)*GSF(1,I)/R + -A(NE+1,NE,4,4)*SF(J)*GSF(2,I)/R END IF ST(IQ,IP) = ST(IQ,IP) + CON3*CNST TANK(IQ,IP) = TANK(IQ,IP)+CON3*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1)+CON4*CNST ST(IQ+1,IP) = ST(IQ+1,IP)+CON4I*CNST TANK(IQ,IP+1) = TANK(IQ,IP+1)+CON4*CNST TANK(IQ+1,IP) = TANK(IQ+1,IP)+CON4I*CNST END IF C------------------------------------------------------------------- C CALCULATE [K33] -- Shear terms C-------------------------------------------------------------------- IQ = I1+2*NIF IP = IN+2*NIF CON5 = A(NE,NE,5,5)*GSF(1,J)*GSF(1,I) + + A(NE,NE,4,5)*(GSF(1,J)*GSF(2,I)+GSF(2,J)*GSF(1,I)) + + A(NE,NE,4,4)*GSF(2,I)*GSF(2,J) IF(IN.LT.JN) THEN CON6 = A(NE,NE+1,5,5)*GSF(1,J)*GSF(1,I) + +A(NE,NE+1,4,5)*(GSF(1,J)*GSF(2,I)+GSF(2,J)*GSF(1,I)) + + A(NE,NE+1,4,4)*GSF(2,I)*GSF(2,J) END IF ST(IQ,IP) = ST(IQ,IP) + CON5*CNST TANK(IQ,IP) = TANK(IQ,IP) + CON5*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1) + CON6*CNST ST(IQ+1,IP) = ST(IQ+1,IP) + CON6*CNST TANK(IQ,IP+1) = TANK(IQ,IP+1) + CON6*CNST TANK(IQ+1,IP) = TANK(IQ+1,IP) + CON6*CNST END IF I1= I1+1 655 CONTINUE 860 JJ = NDF*J+1 770 II = NDF*I+1 600 CONTINUE C---------------------------------------------------------------------- C CARRY OUT REDUCED INTEGRATION FOR TRANSVERSE NORMAL TERMS C (Terms containing Qi3's ,i = 1,2,6) C C---------------------------------------------------------------------- DO 800 NI =1,MGP DO 800 NJ =1,MGP XI = GAUSS(NI,MGP) ETA= GAUSS(NJ,MGP) CALL SHAPE2(NPE,XI,ETA,SF,GSF,DET,ELXY) CNST = DET*WT(NI,MGP)*WT(NJ,MGP) II = 1 DO 870 I = 1,NPE JJ=1 DO 760 J = 1,NPE I1=II JN=JJ+NPLY NE = 0 DO 755 IN = JJ,JN NE = NE+1 C------------------------------------------------------------------- C CALCULATE [K13] -- Transverse normal terms C-------------------------------------------------------------------- IP = IN+2*NIF CON1 = ABAR(NE,NE,1,3)*GSF(1,I)*SF(J) + + ABAR(NE,NE,3,6)*GSF(2,I)*SF(J) IF(IN.LT.JN) THEN CON2 = ABAR(NE,NE+1,1,3)*GSF(1,I)*SF(J) + + ABAR(NE,NE+1,3,6)*GSF(2,I)*SF(J) CON2I= ABAR(NE+1,NE,1,3)*GSF(1,I)*SF(J) + + ABAR(NE+1,NE,3,6)*GSF(2,I)*SF(J) END IF ST(I1,IP) = ST(I1,IP)+CON1*CNST TANK(I1,IP) = TANK(I1,IP)+CON1*CNST IF(IN.LT.JN) THEN ST(I1,IP+1) = ST(I1,IP+1)+CON2*CNST ST(I1+1,IP) = ST(I1+1,IP)+CON2I*CNST TANK(I1,IP+1)= TANK(I1,IP+1)+CON2*CNST TANK(I1+1,IP)= TANK(I1+1,IP)+CON2I*CNST END IF C------------------------------------------------------------------- C CALCULATE [K31] -- Transverse normal terms C-------------------------------------------------------------------- C C ....LINEAR TERMS FOLLOW C IP = IN IQ = I1+2*NIF CON1 = ABAR(NE,NE,1,3)*GSF(1,J)*SF(I) + +ABAR(NE,NE,3,6)*GSF(2,J)*SF(I) IF(IN.LT.JN) THEN CON2 = ABAR(NE+1,NE,1,3)*GSF(1,J)*SF(I) + + ABAR(NE+1,NE,3,6)*GSF(2,J)*SF(I) CON2I= ABAR(NE,NE+1,1,3)*GSF(1,J)*SF(I) + + ABAR(NE,NE+1,3,6)*GSF(2,J)*SF(I) END IF ST(IQ,IP) = ST(IQ,IP) + CON1*CNST TANK(IQ,IP) = TANK(IQ,IP) + CON1*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1) + CON2*CNST ST(IQ+1,IP) = ST(IQ+1,IP) + CON2I*CNST TANK(IQ,IP+1) = TANK(IQ,IP+1) + CON2*CNST TANK(IQ+1,IP) = TANK(IQ+1,IP) + CON2I*CNST END IF C------------------------------------------------------------------- C CALCULATE [K23] -- Transverse normal terms C-------------------------------------------------------------------- IP = IN+2*NIF IQ = I1+NIF CON3 = ABAR(NE,NE,2,3)*GSF(2,I)*SF(J) + + ABAR(NE,NE,3,6)*GSF(1,I)*SF(J) IF(IN.LT.JN) THEN CON4 = ABAR(NE,NE+1,2,3)*GSF(2,I)*SF(J) + + ABAR(NE,NE+1,3,6)*GSF(1,I)*SF(J) CON4I= ABAR(NE+1,NE,2,3)*GSF(2,I)*SF(J) + + ABAR(NE+1,NE,3,6)*GSF(1,I)*SF(J) END IF ST(IQ,IP) = ST(IQ,IP) + CON3*CNST TANK(IQ,IP) = TANK(IQ,IP) + CON3*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1) + CON4*CNST ST(IQ+1,IP) = ST(IQ+1,IP) + CON4I*CNST TANK(IQ,IP+1)= TANK(IQ,IP+1) + CON4*CNST TANK(IQ+1,IP)= TANK(IQ+1,IP) + CON4I*CNST END IF C------------------------------------------------------------------- C CALCULATE [K32] -- Transverse normal terms C-------------------------------------------------------------------- IQ = I1+2*NIF IP = IN+NIF CON3 = ABAR(NE,NE,2,3)*GSF(2,J)*SF(I) + + ABAR(NE,NE,3,6)*GSF(1,J)*SF(I) IF(IN.LT.JN) THEN CON4 = ABAR(NE+1,NE,2,3)*GSF(2,J)*SF(I) + + ABAR(NE+1,NE,3,6)*GSF(1,J)*SF(I) CON4I= ABAR(NE,NE+1,2,3)*GSF(2,J)*SF(I) + + ABAR(NE,NE+1,3,6)*GSF(1,J)*SF(I) END IF ST(IQ,IP) = ST(IQ,IP)+ CON3*CNST TANK(IQ,IP) = TANK(IQ,IP) + CON3*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1) + CON4*CNST ST(IQ+1,IP) = ST(IQ+1,IP) + CON4I*CNST TANK(IQ,IP+1) = TANK(IQ,IP+1) + CON4*CNST TANK(IQ+1,IP) = TANK(IQ+1,IP) + CON4I*CNST END IF C------------------------------------------------------------------- C CALCULATE [K33] -- Transverse normalterms C-------------------------------------------------------------------- IQ = I1+2*NIF IP = IN+2*NIF CON5 = 2*ABAR(NE,NE,2,3)*SF(J)*SF(I)/R + + DBAR(NE,NE,3,3)*SF(J)*SF(I) IF(IN.LT.JN) THEN CON6 = ABAR(NE,NE+1,2,3)*SF(J)*SF(I)/R + + ABAR(NE+1,NE,2,3)*SF(J)*SF(I)/R + + DBAR(NE,NE+1,3,3)*SF(J)*SF(I) END IF ST(IQ,IP) = ST(IQ,IP) + CON5*CNST TANK(IQ,IP) = TANK(IQ,IP) + CON5*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1) + CON6*CNST ST(IQ+1,IP) = ST(IQ+1,IP) + CON6*CNST TANK(IQ,IP+1) = TANK(IQ,IP+1) + CON6*CNST TANK(IQ+1,IP) = TANK(IQ+1,IP) + CON6*CNST END IF I1= I1+1 755 CONTINUE 760 JJ = NDF*J+1 870 II = NDF*I+1 800 CONTINUE SUM =0.0 DO 260 I=1,NN DO 270 J=1,NN SUM = SUM+ST(I,J)*ELFC(J) 270 CONTINUE 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) IF(ISMEAR.EQ.1) THEN CALL REORDR(NN,NPE,NDF,ST,SM,TANK,ELP) END IF RETURN END C_______________________________________________________________________ C C THIS SUBROUTINE CALCULATES THE PLY STIFNESS MATRICES C [A], [ABAR], [DBAR], [D], [DB], [F] AND [G] C C_______________________________________________________________________ SUBROUTINE STFPLY(NPLY,DEP,Z,THICK,THETA,IPLY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NMX=10,NEX=10,NUX=486,NGX=150,NBX=126) COMMON /ABDF/ A(10,10,6,6), ABAR(10,10,6,6),DBAR(10,10,6,6), + G(10,10) COMMON /DDBF/ D(10,10,10,6,6),DB(10,10,10,6,6),F(10,10,10,10,6,6) COMMON /SS/ SQ(15,6,6) DIMENSION Z(15), THICK(15), THETA(20) NIF=NPLY+1 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,10 DO 40 J=1,10 G(I,J) = 0.0 DO 40 K=1,6 DO 40 L=1,6 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,10 D(I,J,M,K,L) = 0.0D0 DB(I,J,M,K,L) = 0.0D0 DO 40 N = 1,10 F(I,J,M,N,K,L)= 0.0D0 40 CONTINUE DO 43 L=1,15 DO 43 I=1,6 DO 43 J=1,6 43 SQ(L,I,J)=0.0 IF(IPLY.EQ.0) CALL CIJ(NPLY) IF(IPLY.EQ.1) CALL CIJS(NPLY,THETA) THK = THICK(1) G(1,1) = THK/3.0D0 DO 50 I =1,6 DO 50 J =1,6 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,6 DO 70 J = 1,6 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,6 DO 90 J = 1,6 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 re-arranges the element stiffness matrices when C the full model for stiffeners is used. (i.e reserves a location C for the stiffness contribution from the stiffeners.) C______________________________________________________________________ SUBROUTINE REORDR(NN,NPE,NDF,ST,SM,TANK,ELP) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER (NMX=10,NEX=10,NUX=486,NGX=150,NBX=126) COMMON /LOCATE/ ILOCS(30), ILOCR(30),ISCOD COMMON /PAR1/ NEIGEN,IFLAG,NPLY,NPLS,NPLR DIMENSION ST(NGX,NGX),SM(NGX,NGX),TANK(NGX,NGX),SS(NGX,NGX) DIMENSION ELP(NGX),ETP(NGX) DO 5 I = 1,NGX 5 ETP(I) = ELP(I) DO 15 I = 1,NGX 15 ELP(I) = 0.0 NIS = NPLY + NPLS + NPLR + 1 NDF = NIS*3 NN = NDF*NPE DO 10 I = 1,NGX DO 10 J = 1,NGX SS(I,J) = 0.0D0 10 ST(I,J) = 0.0D0 NI = NPLR IF(ISCOD.EQ.0) NI = 0 DO 20 I = 1,NPE DO 25 J = 1,3 DO 30 K = 1,NPLY+1 II = (I-1)*NDF + (J-1)*NIS + K + NI IK = (I-1)*(NPLY+1)*3 + (J-1)*(NPLY+1) + K ELP(II) = ETP(IK) DO 35 L = 1,NPE DO 35 M = 1,3 DO 35 N = 1,NPLY+1 JJ = (L-1)*NDF + (M-1)*NIS + N + NI JK = (L-1)*(NPLY+1)*3 + (M-1)*(NPLY+1) + N SS (II,JJ) = SM(IK,JK) ST (II,JJ) = TANK(IK,JK) 35 CONTINUE 30 CONTINUE 25 CONTINUE 20 CONTINUE DO 40 I = 1,NGX DO 40 J = 1,NGX 40 SM(I,J) = SS(I,J) RETURN END C_______________________________________________________________________ C C THIS SUBROUTINE CALCULATES THE PLY STIFFNESS COEFFICIENTS C (Cij's) OF THE SHELL C IN THE GLOBAL DIRECTION C C_______________________________________________________________________ SUBROUTINE CIJ (NPLY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON / SS / SQ(15,6,6) COMMON /CCC/ E1,E2,E3,G12,G13,G23,ANU12,ANU13,ANU23 COMMON /ELAST/ THETA(20),RADIUS(15),ELX(3) PI= DACOS(-1.0D0) ANU21 = ANU12*E2/E1 ANU31 = ANU13*E3/E1 ANU32 = ANU23*E3/E2 DELTA = 1-ANU12*ANU21-ANU13*ANU31-ANU23*ANU32-ANU21*ANU32*ANU13 + -ANU12*ANU23*ANU31 S11 = 1.0/E1 S12 = -ANU21/E2 S13 = -ANU31/E3 S22 = 1.0/E2 S23 = -ANU32/E3 S33 = 1.0/E3 CAPS = S11*S22*S33-S11*S23*S23-S22*S13*S13-S33*S12*S12+ + 2.0*S12*S13*S23 C11 = E1*(1-ANU23*ANU32)/DELTA C12 = E2*(ANU12+ANU13*ANU31)/DELTA C12 = (S13*S23-S12*S33)/CAPS 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 C C C.....CALCULATE THE REDUCED TRANSFORMED STIFFNESS MATRIX C C 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 SQ(L,1,1) = C11*CN4 + 2*(C12+2*C66)*SN2*CN2 +C22*SN4 SQ(L,1,2) = (C11+C22 - 4*C66)* SN2*CN2+C12*(SN4+CN4) SQ(L,1,6) = (C11- C12-2*C66) * SN*CN3 + (C12-C22+2*C66)*SN3*CN SQ(L,2,2) = C11*SN4 + 2*(C12+2*C66)*SN2*CN2 + C22*CN4 SQ(L,2,6) = (C11-C12-2*C66)*SN3*CN + (C12-C22+2*C66)*SN*CN3 SQ(L,6,6) = (C11+C22 -2*C12-2*C66)*SN2*CN2 +C66*(SN4+CN4) SQ(L,2,1) = SQ(L,1,2) SQ(L,6,1) = SQ(L,1,6) SQ(L,6,2) = SQ(L,2,6) SQ(L,1,3) = C13*CN2 + C23*SN2 SQ(L,2,3) = C13*SN2 + C23*CN2 SQ(L,3,3) = C33 SQ(L,3,6) = (C13-C23)*CN*SN SQ(L,4,4) = G23*CN2 + G13*SN2 SQ(L,4,5) = (G13-G23) *CN*SN SQ(L,5,5) = G13*CN2 +G23*SN2 SQ(L,3,1) = SQ(L,1,3) SQ(L,3,2) = SQ(L,2,3) SQ(L,6,3) = SQ(L,3,6) SQ(L,5,4) = SQ(L,4,5) 30 CONTINUE RETURN END C---------------------------------------------------------------------- C C This Subroutine calculates the ply stiffness coefficients of C the stiffeners in the global direction. C----------------------------------------------------------------------- SUBROUTINE CIJS (NPLY,THETA) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON / SS / SQ(15,6,6) COMMON /CCC/ E1,E2,E3,G12,G13,G23,ANU12,ANU13,ANU23 DIMENSION THETA(20) PI = DACOS(-1.0D0) 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,4,4) = (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 C THIS SUBROUTINE CALCULATES THE LAYERWISE ELEMENT DIRECT AND C TANGENT STIFFNESS MATRIX OF A STIFFENER ELEMENT OF A DISCRETELY C STIFFENED LAMINATED SHELL OR PLATE. C C_______________________________________________________________________ SUBROUTINE STFFUL(NPE,RHO,NDB,R,NRED,IPR,N,ILOC,THT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NMX=10,NEX=10,NUX=486,NGX=150,NBX=126) COMMON/ABDF/ A(10,10,6,6), ABAR(10,10,6,6), DBAR(10,10,6,6), + G(10,10) COMMON/DDBF/ D(10,10,10,6,6),DB(10,10,10,6,6),F(10,10,10,10,6,6) COMMON/DOUBLE/ DVDS(20),DWDS(20),TANK(NGX,NGX),DUDX(20),DUDY(20) COMMON/GLOBAL/ GSTIF(NUX,NBX),GF(NUX),W(NUX),EPL(NUX) COMMON/SMAS/ GMASS(NEX,NEX),EGNVAL(NEX) COMMON/PAR1/ NEIGEN,IFLAG,NPLY,NPLS,NPLR COMMON/SEL/ ST(NGX,NGX),ELP(NGX),ELXY(9,2),ELFC(NGX) COMMON/DEPTH/ Z(15),ZZ(15),THICK(15),SM(NGX,NGX) COMMON/SHP/ SF(4),GSF(4),GDDSF(4),DET COMMON/ELAST/ THETA(20),RADIUS(15),ELX(3) COMMON/MSR/ NODR(50,4),NODS(50,4) COMMON/DIMF/ DWKDS(20) DIMENSION WT(4,4),GAUSS(4,4),ILOC(10) 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 IPROB = IPR K = 1 NIF = NPLS+1 L = 4 NIS = NPLY+1 NIT = NPLY + NPLS + NPLR + 1 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.4) NGP = NGP-1 DO 5 I = 1,NGX DO 5 J = 1,NGX TANK(I,J) = 0.0D0 SM(I,J) = 0.0D0 5 ST(I,J) = 0.0D0 DO 10 I = 1,NGX 10 ELP(I) = 0.0D0 H1 = DABS(ELX(NPE)-ELX(1)) C_______________________________________________________________________ C C......FULL INTEGRATION OF TANGENT STIFFNESS MATRIX FOLLOWS C C_______________________________________________________________________ DO 100 NI = 1,NGP XI = GAUSS(NI,NGP) CALL SHAPE1(XI,H1,NPE,NET,IELEM) CNST = DET*WT(NI,NGP) DO 15 I = 1,20 DVDS(I) = 0.0 DWDS(I) = 0.0 15 DWKDS(I) = 0.0 DO 130 I=1,NPE L1 =(I-1)*NDF+1+NIF L2 =(I-1)*NDF+1 L3 =(I-1)*NDF+1+2*NIF DO 140 J= 1,NIF DVDS(J) = DVDS(J)+GSF(I)*ELFC(L2) DWDS(J) = DWDS(J)+GSF(I)*ELFC(L1) DWKDS(J) =DWKDS(J)+GSF(I)*ELFC(L3) L1 = L1+1 L2 = L2+1 L3 = L3+1 140 CONTINUE 130 CONTINUE II=1 DO 70 I=1,NPE JJ=1 DO 60 J=1,NPE I1=II JN=JJ+NPLS 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 CON1D = 0.0 CON2D = 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 CON55M = 0.0 CON66M = 0.0 CON55S = 0.0 CON66S = 0.0 C_______________________________________________________________________ C C..... CALCULATE K11(TANGENT AND DIRECT) C C_______________________________________________________________________ ST(I1,IN) = ST(I1,IN) +(A(NE,NE,K,K)*GSF(I)*GSF(J) + + DBAR(NE,NE,L,L)*SF(I)*SF(J) + - ABAR(NE,NE,L,L)*SF(I)*SF(J)/R + - ABAR(NE,NE,L,L)*SF(I)*SF(J)/R + + A(NE,NE,L,L)*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,K,K)*GSF(I)*GSF(J) + +DBAR(NE,NE+1,L,L)*SF(I)*SF(J) + -ABAR(NE,NE+1,L,L)*SF(I)*SF(J)/R + -ABAR(NE+1,NE,L,L)*SF(I)*SF(J)/R + +A(NE,NE+1,L,L)*SF(I)*SF(J)/R/R)*CNST ST(I1+1,IN) = ST(I1+1,IN) +(A(NE+1,NE,K,K)*GSF(I)*GSF(J) + +DBAR(NE+1,NE,L,L)*SF(I)*SF(J) + -ABAR(NE,NE+1,L,L)*SF(I)*SF(J)/R + -ABAR(NE+1,NE,L,L)*SF(I)*SF(J)/R + +A(NE+1,NE,L,L)*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..... CALCULATE K12 (TANGENT AND DIRECT) C C_______________________________________________________________________ C ....LINEAR TERMS FOLLOW IP = IN+NIF CON1 = A(NE,NE,K,K)*GSF(I)*SF(J)/R + + ABAR(NE,NE,L,L)*SF(I)*GSF(J) - + A(NE,NE,L,L)*SF(I)*GSF(J)/R + + ABAR(NE,NE,1,3)*SF(J)*GSF(I) IF(IN.LT.JN) THEN CON2 = A(NE,NE+1,K,K)*GSF(I)*SF(J)/R + + ABAR(NE+1,NE,L,L)*SF(I)*GSF(J) - + A(NE,NE+1,L,L)*SF(I)*GSF(J)/R + + ABAR(NE,NE+1,1,3)*SF(J)*GSF(I) CON2I= A(NE+1,NE,K,K)*GSF(I)*SF(J)/R + + ABAR(NE,NE+1,L,L)*SF(I)*GSF(J) - + A(NE+1,NE,L,L)*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,K,K)*GSF(I)*GSF(J)*(DWDS(ME)+DWKDS(ME)) CON1D=CON1D+D(NE,NE,ME,K,K)*GSF(I)*GSF(J)*(DWDS(ME)+2*DWKDS(ME))/2 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,K,K)*GSF(I)*GSF(J)*(DWDS(IZ)+DWKDS(IZ)) CON2D=CON2D+D(NE,NE+1,IZ,K,K)*GSF(I)*GSF(J)*(DWDS(IZ)+ + 2*DWKDS(IZ))/2 IZ = IZ+1 160 CONTINUE END IF ST(I1,IP) = ST(I1,IP)+(CON1+CON1D)*CNST TANK(I1,IP) = TANK(I1,IP)+(CON1+CON11)*CNST IF(IN.LT.JN) THEN ST(I1,IP+1) = ST(I1,IP+1)+(CON2+CON2D)*CNST ST(I1+1,IP) = ST(I1+1,IP)+(CON2I+CON2D)*CNST TANK(I1,IP+1) =TANK(I1,IP+1)+(CON2+CON22)*CNST TANK(I1+1,IP) =TANK(I1+1,IP)+(CON2I+CON22)*CNST END IF C_______________________________________________________________________ C C CALCULATE K21 (TANGENT AND DIRECT) C C_______________________________________________________________________ IP = IN IQ = I1+NIF CON1 = A(NE,NE,K,K)*GSF(J)*SF(I)/R + + ABAR(NE,NE,L,L)*SF(J)*GSF(I) - + A(NE,NE,L,L)*SF(J)*GSF(I)/R + + ABAR(NE,NE,1,3)*SF(I)*GSF(J) IF(IN.LT.JN) THEN CON2 = A(NE,NE+1,K,K)*GSF(J)*SF(I)/R + + ABAR(NE,NE+1,L,L)*SF(J)*GSF(I) - + A(NE,NE+1,L,L)*SF(J)*GSF(I)/R + + ABAR(NE+1,NE,1,3)*SF(I)*GSF(J) CON2I= A(NE+1,NE,K,K)*GSF(J)*SF(I)/R + + ABAR(NE+1,NE,L,L)*SF(J)*GSF(I) - + A(NE,NE+1,L,L)*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,K,K)*GSF(J)*GSF(I)*(DWDS(ME)+DWKDS(ME)) 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,K,K)*GSF(J)*GSF(I)*(DWDS(ME)+DWKDS(ME)) IZ = IZ+1 196 CONTINUE END IF ST(IQ,IP) = ST(IQ,IP)+(CON1+CON33)*CNST TANK(IQ,IP) = TANK(IQ,IP)+(CON1+CON33)*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1)=ST(IQ,IP+1)+(CON2+CON44)*CNST ST(IQ+1,IP)=ST(IQ+1,IP)+(CON2I+CON44)*CNST TANK(IQ,IP+1) =TANK(IQ,IP+1)+(CON2+CON44)*CNST TANK(IQ+1,IP) =TANK(IQ+1,IP)+(CON2I+CON44)*CNST END IF C_______________________________________________________________________ C C CALCULATE K22 (TANGENT AND DIRECT) C C_______________________________________________________________________ C ....LINEAR TERMS FOLLOW IQ = I1+NIF IP = IN+NIF CON5 = A(NE,NE,K,K)*SF(I)*SF(J)/R/R + + A(NE,NE,L,L)*GSF(I)*GSF(J) + +ABAR(NE,NE,1,3)*SF(I)*SF(J)/R + +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,K,K)*SF(I)*SF(J)/R/R +