PROGRAM WIRE3
C======================================================================
C----- UNEVEN ELEMENT<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C---------DZMAXER = MAX ERROR IN Z DURING ITERATIVE PROCEDURE
C DZMAXER IS DIVIDED BY DOMAIN LENGTH
C TX = TENSION IN X-DIRECTION
C TX IS DIVIDED BY GRAVITY.
C======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=2, MXE=1000, MXN=MXE+1,NBW=2*(ND-1)+1 )
PARAMETER ( MXITERA=80, ERRMAX=1.D-4, GRAVITY = 1.D0 )
DIMENSION NODEX(MXE,ND),X(MXN),Z(MXN),PARA(MXN),
* STIFF(ND,ND),A(MXN,NBW),RHS(MXN),TS(MXN)
C======================================================================
C (1)-READING OF DATA
CALL INPUT (NE,DOMAINL,WDENSITY,TXMAX,TXMIN,TXSTEP, Z1,ZN)
GW = WDENSITY*GRAVITY
C======================================================================
C (2)-FINITE ELEMENT CREATION
CALL ELEMENT ( DOMAINL, MXE,MXN,ND,NE,NNODE,NODEX,X )
C======================================================================
C (4)-COMPUTATION OF INITIAL LOCATION OF STRING
C-------INITIAL Z STANDS FOR A SET OF GUESS STRING LOCATION.
CALL INITIALZ ( MXN,NNODE,GW,TXMAX,X, Z1, ZN, Z )
C======================================================================
OPEN ( 1,FILE='SOLUTION.FEM', STATUS='UNKNOWN')
OPEN ( 2,FILE='SUMMARY.FEM', STATUS='UNKNOWN')
WRITE(1,*) 'X Z PARABOLIC TENSION TX ITERATIONS'
WRITE(2,*) 'TX TS-UP TS-DW MAXDEPTH L-OF-WIRE W-OF-WIRE SUM-TZ',
* ' %ERROR ITERATIONS'
C
C
NSTEPS = (TXMAX-TXMIN)/TXSTEP + 1
TX = TXMAX + TXSTEP
DO NS = 1, NSTEPS
TX = TX - TXSTEP
IF ( TX .LE. 0.D0 ) EXIT
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
ITERA = 0
DZMAXER = 2.D0*ERRMAX
DO WHILE ( ITERA .LE. MXITERA .AND. DZMAXER .GT. ERRMAX )
C*********************************************************************
ITERA = ITERA + 1
C==============================================
C (7)-CONSTRUCTION OF FEM-MATRIX EQUATION
CALL MATRIX ( MXE,MXN,ND,NBW,NE,NNODE, GRAVITY,WDENSITY,
* TX,STIFF, NODEX, X, Z, A,RHS )
C==============================================
C (8)-IMPLEMENTATION OF BOUNDARY CONDITION
CALL FORM ( MXN, NBW, NNODE, A, RHS, Z1, ZN )
C==============================================
C (9)-SOLVE FOR UNKNOWN VARIABLES
CALL SYSTEM ( MXN, NBW, NNODE, A , RHS )
C==============================================
C(10)-EVALUATION OF NEW STRING LOCATION BASED UPON OLD U AND NEW U.
C-------- UMAX = MAXIMUM ERROR IN U(I) DURING ITERATION
CALL NEWZ ( MXN,NNODE,DOMAINL, RHS ,Z, DZMAXER )
END DO
C-------- END OF ITARETIVE PROCEDURE ---------
C*********************************************************************
CALL PRT (ND,MXE,MXN,NE,NNODE,NODEX,X,TX,GW,PARA,RHS,TS,
* DZMAXER,ITERA)
C------ END OF COMPUTATIONAL PROCEDURE
END DO
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
CLOSE (2)
CLOSE (1)
STOP'NORMAL TERMINATION'
END
C
C
C23456789012345678901234567890123456789012345678901234567890123456789012
SUBROUTINE PRT (ND,MXE,MXN,NE,NNODE,NODEX,X,TX,GW,PARA,RHS,TS,
* DZMAXER,ITERA)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION X(MXN),Z(MXN),PARA(MXN),RHS(MXN),TS(MXN),NODEX(MXE,ND)
CALL PARAB (MXN,NNODE,X,TX,GW,RHS,PARA)
CALL CMPDEPTH ( MXN,NNODE,X,RHS,DEPTH,WL )
CALL TENSION (ND,MXE,MXN,NE,NODEX,X,RHS,TX,TS,SUMTZ)
WOW = GW*WL
WRITE(1,*)
WRITE(1,*) X(1),RHS(1), PARA(1), TS(1),TX,ITERA
DO I = 2 , NNODE
WRITE(1,*) X(I),RHS(I), PARA(I), TS(I)
END DO
PRCNTER = ( WOW - SUMTZ ) / WOW * 100.D0
WRITE (2,*) TX,TS(1),TS(NNODE), DEPTH, WL, WOW,SUMTZ,PRCNTER,
* ITERA
RETURN
END
C
C
SUBROUTINE PARAB (MXN,NNODE,X,TX,GW,RHS,PARA)
C------- COMPUTES RADIUS AT TX INFINITE STRENGTH.-------
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION X(MXN), PARA(MXN),RHS(MXN)
CONST = 0.5D0*GW/TX
TL = X(NNODE) - X(1)
C = (1.D0/TL)*(RHS(NNODE)-CONST*TL*TL)
DO I = 1 , NNODE
PARA(I) = CONST*X(I)*X(I)+C*X(I)+RHS(1)
END DO
RETURN
END
C
C
SUBROUTINE CMPDEPTH (MXN,NNODE,X,RHS,DEPTH,WL )
C------- COMPUTES MAX. DEPTH AND LENGTH OF WIRE -------
C------- DEPTH = DEPTH OF WIRE FROM STRAIGHT LINE CONNECTING BOUNDARY.
C------- WL = LENGTH OF WIRE
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION X(MXN), RHS(MXN)
WL = 0.D0
SLOPE = (RHS(NNODE) - RHS(1))/(X(NNODE)-X(1))
DEPTH = 0.D0
DO I = 2 , NNODE
Z0 = SLOPE*(X(I)-X(1)) + RHS(1)
ZX = RHS(I) - Z0
DEPTH = DMAX1 ( DEPTH, DABS(ZX) )
DX = X(I) - X(I-1)
DZ = RHS(I) - RHS(I-1)
WL = WL + DSQRT ( DX*DX + DZ*DZ )
END DO
RETURN
END
C
C
SUBROUTINE TENSION (ND,MXE,MXN,NE,NODEX,X,RHS,TX,TS,SUMTZ)
C------- COMPUTES RADIUS AT TX INFINITE STRENGTH.-------
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION X(MXN), RHS(MXN), TS(MXN), NODEX(MXE,ND)
C-------- TENSION IN EACH ELEMENT
NNODE = NE + 1
DO I = 1 , NNODE
TS(I) = 0.D0
END DO
DO IEL = 1 , NE
I = NODEX(IEL,1)
J = NODEX(IEL,2)
DX = X(J) - X(I)
DZ = DABS (RHS(J) - RHS(I))
DS = DSQRT( DX**2 + DZ**2 )
TS(I) = TS(I) + TX*DS/DX
TS(J) = TS(J) + TX*DS/DX
END DO
DO I = 2 , NE
TS(I) = TS(I) / 2.D0
END DO
C-------- UPSTREAM
DX = X(2) - X(1)
DZ = DABS ( RHS(2) - RHS(1) )
SUMTZ = TX*DZ/DX
C-------- DWONSTREAM
DX = X(NNODE) - X(NNODE-1)
DZ = DABS ( RHS(NNODE) - RHS(NNODE-1) )
SUMTZ = SUMTZ + TX*DZ/DX
RETURN
END
C
C
SUBROUTINE MATRIX ( MXE,MXN,ND,NBW,NE,NNODE, GRAVITY,WDENSITY,
* TX,STIFF, NODEX, X, Z, A,RHS )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),X(MXN),Z(MXN),RHS(MXN),
* STIFF(ND,ND),
* A(MXN,NBW)
CONST = GRAVITY * WDENSITY / TX
DO I = 1 , NNODE
A(I,1) = 0.D0
A(I,2) = 0.D0
A(I,3) = 0.D0
RHS(I) = 0.D0
END DO
DO IEL = 1 , NE
I = NODEX(IEL,1)
J = NODEX(IEL,2)
DX = X(J) - X(I)
DZ = Z(J) - Z(I)
DSDX = DSQRT ( 1.D0 + (DZ/DX)**2 )
FORCE = CONST * DSDX * DX / 2.D0
EE = 1.D0/DX
STIFF(1,1) = -EE
STIFF(1,2) = EE
STIFF(2,1) = EE
STIFF(2,2) = -EE
A(I,2) = A(I,2) + STIFF(1,1)
A(I,3) = A(I,3) + STIFF(1,2)
A(J,1) = A(J,1) + STIFF(2,1)
A(J,2) = A(J,2) + STIFF(2,2)
RHS(I) = RHS(I) + FORCE
RHS(J) = RHS(J) + FORCE
END DO
RETURN
END
C
C
SUBROUTINE INITIALZ ( MXN,NNODE,GW,TXMAX,X, Z1, ZN, Z )
C------- COMPUTES RADIUS AT TX INFINITE STRENGTH.-------
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION X(MXN), Z(MXN)
CONST = 0.5D0*GW/TXMAX
TL = X(NNODE) - X(1)
C = (1.D0/TL)*(ZN-CONST*TL*TL)
DO I = 1 , NNODE
Z(I) = CONST*X(I)*X(I)+C*X(I)+Z1
END DO
RETURN
END
C
C
SUBROUTINE NEWZ ( MXN,NNODE,DOMAINL, RHS ,Z, DZMAXER )
C------- COMPUTES RADIUS AT TX INFINITE STRENGTH.-------
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION RHS(MXN), Z(MXN)
DZMAXER = DABS ( RHS(1) - Z(1) ) / DOMAINL
DO I = 1 , NNODE
DZMAXER = DMAX1 ( DZMAXER, DABS(RHS(I)-Z(I))/DOMAINL )
Z(I) = (RHS(I) + Z(I))/2.D0
END DO
RETURN
END
C
C
SUBROUTINE INPUT (NE,DOMAINL,WDENSITY,TXMAX,TXMIN,TXSTEP, Z1,ZN)
IMPLICIT REAL*8 ( A-H , O-Z )
COMPLEX*16 UPSTRM, DWSTRM
OPEN ( 1,FILE='WIRE.DAT', STATUS='OLD')
C-------- DOMAINL = LENGTH OF DOMAIN IN Z
C-------- WDENSITY = DENSITY OF WIRE AS GRAM PER CENTIMETER
C-------- TXMAX = POSSIBLE MAX. TENSION IN X-DIRECTION
C-------- YXMIN = POSSIBLE MIN. TENSION IN X-DIRECTION
C-------- TXSTEP = INCREMENTAL TENSION
C-------- Z1 = ELEVATION OF WIRE AT X=0
C-------- ZN = ELEVATION OF WIRE AT X=DOMAINL
C------- READING CONTROL PARAMETERS ------------------
READ(1,*) NE
READ(1,*) DOMAINL
READ(1,*) WDENSITY
READ(1,*) TXMAX,TXMIN,TXSTEP
READ(1,*) Z1
READ(1,*) ZN
CLOSE (1)
RETURN
END
C
C
SUBROUTINE ELEMENT ( DOMAINL,MXE,MXN,ND,NE,NNODE,NODEX,X )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),X(MXN)
C--------NE = NUMBER OF ELEMENT
SLOPE = 0.D0
A = 2.D0*(SLOPE-1.D0)
B = 3.D0*(1.D0-SLOPE)
DO I = 1 , NE
NODEX(I,1) = I
NODEX(I,2) = I+1
END DO
NNODE = NE + 1
DX = DOMAINL / NE
X(1) = 0.
DO I = 2 , NNODE
T = (I-1)*DX/DOMAINL
X(I) = DOMAINL * T*(T*(A*T + B) + SLOPE)
END DO
RETURN
END
C
C
SUBROUTINE FORM ( MXN, NBW, NNODE, A, RHS, Z1, ZN )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A(MXN,NBW), RHS(MXN)
C-------UP STREAM SIDE
A(1,2) = 1.D0
A(1,3) = 0.D0
RHS(1) = Z1
RHS(2) = RHS(2) - Z1*A(2,1)
A(2,1) = 0.D0
C-------DOWN STREAM SIDE
A(NNODE,2) = 1.D0
A(NNODE,1) = 0.D0
RHS(NNODE) = ZN
RHS(NNODE-1) = RHS(NNODE-1) - ZN*A(NNODE-1,3)
A(NNODE-1,3) = 0.D0
RETURN
END
C
C
SUBROUTINE SYSTEM ( MXN, NBW, NNODE, A , B )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A(MXN,NBW) , B(MXN)
B(1) = B(1) / A(1,2)
A(1,2) = A(1,3) / A(1,2)
DO I = 2 , NNODE
P = A(I,2) - A(I,1) * A(I-1,2)
A(I,2) = A(I,3) / P
B(I) = ( B(I) - A(I,1)*B(I-1) ) / P
END DO
C------ BACK SUBSTITUTION ----
DO I = NNODE-1, 1,-1
B(I) = B(I) - A(I,2) * B(I+1)
END DO
RETURN
END