PROGRAM POTENTIAL
COMMON
1/GRID/XOR,YOR,DELTAX,DELTAY,X(1500),Y(1500),NXI,NYI,XMIN,YMIN
1 ,XMAX,YMAX,THETA(1500),R(1500,1500)
1/KINDPOT/POT(1500,1500),POT1(1500,1500)
1/INUMFLG/IFLAG1,IFLAG2,IFLAG3,IFLAG4
1/PARFLOW/U0
1/SOURSINK/XORT(4),YORT(4),Q(4),B(4)
1/VEL/U(1500,1500),V(1500,1500),UTOT(1500,1500)
1/PRES/PR(1500,1500)
1/IFL/NNTOTI,NNTOTJ
1/ICOORD/IGRD
CHARACTER TITLE*40,VAR*200,TITLE1*40,VAR1*200
LOGICAL IPOT1,IPOT2,IPOT3,IPOT4
INTEGER IRESTART
C-----OPEN & READ DATA FOR THE PROBLEM
C-----GEOMETRY
OPEN(10,FILE='GEOMETRY.IN')
READ(10,100)XMIN,YMIN,XMAX,YMAX,XOR,YOR
READ(10,110)NXI,NYI
READ(10,110)IRESTART
READ(10,110)IGRD
C-----CHECK FOR RESTART OR NOT
CIRESTART=1 ---RESTART
CIRESTART=2 ---NOT RESTART
IF(IRESTART.EQ.1)CALL RESTART
CIF(IRESTART.EQ.1)CALL RESTART
C-----SET THE KIND OF FLOWS
OPEN(30,FILE='DATA.IN')
READ(30,120)IPOT1,IPOT2,IPOT3,IPOT4
READ(30,110)IFLAG1
READ(30,100)U0
READ(30,110)IFLAG2
READ(30,100)Q(1),B(1),XORT(1),YORT(1),Q(2),B(2),XORT(2),YORT(2)
1 ,Q(3),B(3),XORT(3),YORT(3),Q(4),B(4),XORT(4),YORT(4)
100FORMAT(9X,F12.5)
110FORMAT(9X,I6)
120FORMAT(9X,L8)
C-----INITIAL VALUES OF POTETIONAL
IF(IRESTART.EQ.1)GOTO 500
DO I=1,NXI
DO J=1,NYI
POT(I,J)=0.0
U(I,J)=0.0
V(I,J)=0.0
UTOT(I,J)=0.0
PR(I,J)=0.0
ENDDO
ENDDO
C-----BUILD X & Y GRID LINES
500CALL CARTGRID
C-----SET THE KIND OF FLOWS
IF(IPOT1)CALL POTE1
IF(IPOT2)CALL POTE2
IF(IPOT3)CALL POTE3
IF(IPOT4)CALL POTE4
C-----COMPUTE VELOCITY FIELD
CALL VELCOM
C-----COMPUTE PRESSURE FIELD
CALL PRESCOM
C-----RESTART FILE
CALL RESTART
C-----OUTPUT FILE
999OPEN(20,FILE='POTENTIONAL.PLT')
TITLE=' TITLE=" POTENTIONAL FLOW "'
VAR=' VARIABLES = "X","Y","THETA",'//'"POTENTIONAL","U","V",
1 "UTOT"'
c,"PR"'
WRITE(20,*) TITLE
WRITE(20,*) VAR
WRITE(20,201) NXI,NYI
DO J=1,NYI
DO I=1,NXI
WRITE(20,200) X(I),Y(J),THETA(J),POT(I,J),U(I,J),V(I,J)
1 ,UTOT(I,J)
c,PR(I,J)
ENDDO
ENDDO
CLOSE(20)
WRITE(*,*)'TECPLOT FILE IS SAVED:POTENTIONAL.PLT'
200FORMAT(7(1X,F12.5))
201FORMAT(' ZONE I=',I7,',J=',I5,',K= 1, F=POINT')
END
C-----SUBROUTINE CARTGRID
C *******PURPOSE
C MAKES CARTESIAN GRID LINES
C CDD/2001
SUBROUTINE CARTGRID
COMMON
1/GRID/XOR,YOR,DELTAX,DELTAY,X(1500),Y(1500),NXI,NYI,XMIN,YMIN
1 ,XMAX,YMAX,THETA(1500),R(1500,1500)
1/SOURSINK/XORT(4),YORT(4),Q(4),B(4)
1/INUMFLG/IFLAG1,IFLAG2,IFLAG3,IFLAG4
1/IFL/NNTOTI,NNTOTJ
1/ICOORD/IGRD
IF(IGRD.EQ.1)GOTO 91
IF(IGRD.EQ.2)GOTO 92
C-----X GRID LINES
CX(1)=XOR
91DELTAX=(ABS(XMAX)+ABS(XMIN))/FLOAT((NXI-1))
X(1)=XMIN
CDO K=1,IFLAG2
CXORT(K)=XOR
DO I=2,NXI-1
X(I)=X(I-1)+DELTAX
ENDDO
CENDDO
X(NXI)=XMAX
X(NXI+1)=XMAX+DELTAX
C-----Y GRID LINES
CY(1)=YOR
DELTAY=(ABS(YMAX)+ABS(YMIN))/FLOAT((NYI-1))
Y(1)=YMIN
DO J=2,NYI-1
Y(J)=Y(J-1)+DELTAY
ENDDO
Y(NYI)=YMAX
Y(NXI+1)=YMAX+DELTAY
GOTO 93
C-----POLAR COORDINATES
92PI=4.0*ATAN(1.0)
QQ=(XMAX-XOR)/(NXI-1)
AA=2.0*PI/(NYI-1)
THETA(1)=0
DO J=2,NYI
THETA(J)=THETA(J-1)+AA
ENDDO
c-----x/r
X(1)=XOR
DELTAX=XMAX/FLOAT(NXI-1)
DO I=1,NXI
X(I+1)=X(I)+DELTAX
ENDDO
X(NXI+1)=XMAX+DELTAX
93RETURN
END
SUBROUTINE POTE
C-----SUBROUTINE CARTGRID
C *******PURPOSE
C COMPUTES THE STREAM FUNCTION FOR DIFFERENT TYPES OF FLOW
C CDD/2001
COMMON
1/GRID/XOR,YOR,DELTAX,DELTAY,X(1500),Y(1500),NXI,NYI,XMIN,YMIN
1 ,XMAX,YMAX,THETA(1500),R(1500,1500)
1/KINDPOT/POT(1500,1500),POT1(1500,1500)
1/INUMFLG/IFLAG1,IFLAG2,IFLAG3,IFLAG4
1/PARFLOW/U0
1/SOURSINK/XORT(4),YORT(4),Q(4),B(4)
1/IFL/NNTOTI,NNTOTJ
1/ICOORD/IGRD
REAL AL
C-----POT1-PARALLEL FLOW
ENTRY POTE1
DO I=1,NXI+1
DO J=1,NYI+1
IF(IGRD.EQ.2)THEN
POT1(I,J)=U0*X(I)*SIN(THETA(J))
ELSE
POT1(I,J)=U0*Y(J)
ENDIF
ENDDO
ENDDO
CALL TOTALPOT
RETURN
C-----POT2-SOURCE OR SINK
ENTRY POTE2
DO K=1,IFLAG2
AA=Q(K)/(2*3.14*B(K))
XOR=XORT(K)
YOR=YORT(K)
DO J=1,NYI+1
DO I=1,NXI+1
IF(IGRD.EQ.2)THEN
POT1(I,J)=AA*(theta(j))
ELSE
POT1(I,J)=AA*ATAN((Y(J)-YOR)/(X(I)-XOR))
ENDIF
ENDDO
ENDDO
CALL TOTALPOT
ENDDO
RETURN
C-----POT3-CCW LINE VORTEX
ENTRY POTE3
DO K=1,IFLAG2
AA=Q(K)/(2.*3.14*B(K))
XOR=XORT(K)
YOR=YORT(K)
DO I=1,NXI+1
DO J=1,NYI+1
IF(IGRD.EQ.2)THEN
IF((X(I)-XOR).LT.0.001)THEN
POT1(I,J)=0.0
ELSE
POT1(I,J)=AA*log(X(I))
ENDIF
ELSE
POT1(I,J)=AA*LOG(SQRT((X(I)-XOR)**2+(Y(J)-YOR)**2))
ENDIF
ENDDO
ENDDO
CALL TOTALPOT
ENDDO
RETURN
C-----POT4
ENTRY POTE4
reysta=(1000*9.81*0.866*0.008*0.008)/(2*0.001)
do i=1,nxi+1
do j=1,nyi+1
c------
cpot1(i,j)=753.6*LOG(SQRT(X(I)**2+Y(J)**2))
c 1 +(25.0)*(y(j))*(1-4.0/(x(i)**2+y(j)**2))
c -753.6*ATAN(Y(J)/X(I))
c 1 +(25.0)*(x(i))*(1+4.0/(x(i)**2+y(j)**2))
c-----stream function
c------
c-----reysta 17/2/2002
pot1(i,j)=-reysta*y(j)*(1-(1.0/3.0)*(y(j)/0.008)*(y(j)/0.008))
pot1(i,j)=pot1(i,j)+(1000*9.81*0.866*0.008*0.008*0.008)/(3*0.001)
enddo
enddo
CALL TOTALPOT
RETURN
END
SUBROUTINE TOTALPOT
C-----SUBROUTINE TOTALPOT
C *******PURPOSE
C SUMS THE TOTAL STREAM FUNCTION
C CDD/2001
COMMON
1/GRID/XOR,YOR,DELTAX,DELTAY,X(1500),Y(1500),NXI,NYI,XMIN,YMIN
1 ,XMAX,YMAX,THETA(1500),R(1500,1500)
1/KINDPOT/POT(1500,1500),POT1(1500,1500)
1/INUMFLG/IFLAG1,IFLAG2,IFLAG3,IFLAG4
1/PARFLOW/U0
1/SOURSINK/XORT(4),YORT(4),Q(4),B(4)
1/IFL/NNTOTI,NNTOTJ
1/ICOORD/IGRD
DO I=1,NXI+1
DO J=1,NYI+1
POT(I,J)=POT(I,J)+POT1(I,J)
ENDDO
ENDDO
RETURN
END
SUBROUTINE RESTART
C-----SUBROUTINE RESTART
C *******PURPOSE
C MODIDY THE STREAM LINES
C CDD/2001
COMMON
1/GRID/XOR,YOR,DELTAX,DELTAY,X(1500),Y(1500),NXI,NYI,XMIN,YMIN
1 ,XMAX,YMAX,THETA(1500),R(1500,1500)
1/KINDPOT/POT(1500,1500),POT1(1500,1500)
1/INUMFLG/IFLAG1,IFLAG2,IFLAG3,IFLAG4
1/PARFLOW/U0
1/SOURSINK/XORT(4),YORT(4),Q(4),B(4)
1/VEL/U(1500,1500),V(1500,1500),UTOT(1500,1500)
1/IFL/NNTOTI,NNTOTJ
1/ICOORD/IGRD
C-----
OPEN(88,FILE='RESTART.TXT')
DO J=1,NYI
DO I=1,NXI
WRITE(88,300) X(I),Y(J),POT(I,J)
ENDDO
ENDDO
CLOSE(88)
300FORMAT(3(1X,F12.5))
RETURN
END
SUBROUTINE VELCOM
C-----SUBROUTINE VELCOM
C *******PURPOSE
C COMPUTES THE VELOCITIES
C CDD/2001
COMMON
1/GRID/XOR,YOR,DELTAX,DELTAY,X(1500),Y(1500),NXI,NYI,XMIN,YMIN
1 ,XMAX,YMAX,THETA(1500),R(1500,1500)
1/KINDPOT/POT(1500,1500),POT1(1500,1500)
1/INUMFLG/IFLAG1,IFLAG2,IFLAG3,IFLAG4
1/PARFLOW/U0
1/SOURSINK/XORT(4),YORT(4),Q(4),B(4)
1/VEL/U(1500,1500),V(1500,1500),UTOT(1500,1500)
1/IFL/NNTOTI,NNTOTJ
1/ICOORD/IGRD
C-----U VELOCITY
DO I=1,NXI
DO J=1,NYI
IF(IGRD.EQ.2)THEN
IF(X(I).EQ.0.0)THEN
U(I,J)=0.0
ELSE
U(I,J)=1./X(I)*(POT(I,J+1)-POT(I,J))/(THETA(J+1)-THETA(J))
ENDIF
ELSE
U(I,J)=(POT(I,J+1)-POT(I,J))/(Y(J+1)-Y(J))
ENDIF
ENDDO
ENDDO
C-----V VELOCITY
DO J=1,NYI
DO I=1,NXI
V(I,J)=-(POT(I+1,J)-POT(I,J))/(X(I+1)-X(I))
ENDDO
ENDDO
C-----TOTAL VELOCITY
DO J=1,NYI
DO I=1,NXI
UTOT(I,J)=SQRT(U(I,J)**2+V(I,J)**2)
ENDDO
ENDDO
RETURN
END
SUBROUTINE PRESCOM
C-----SUBROUTINE PRESCOM
C *******PURPOSE
C COMPUTES THE PRESSURE FIELD USING BERNOULI EQUATION
C CDD/2001
COMMON
1/GRID/XOR,YOR,DELTAX,DELTAY,X(1500),Y(1500),NXI,NYI,XMIN,YMIN
1 ,XMAX,YMAX,THETA(1500),R(1500,1500)
1/KINDPOT/POT(1500,1500),POT1(1500,1500)
1/INUMFLG/IFLAG1,IFLAG2,IFLAG3,IFLAG4
1/PARFLOW/U0
1/SOURSINK/XORT(4),YORT(4),Q(4),B(4)
1/VEL/U(1500,1500),V(1500,1500),UTOT(1500,1500)
1/PRES/PR(1500,1500)
1/IFL/NNTOTI,NNTOTJ
1/ICOORD/IGRD
REAL P0,UPR
C-----PRESSURE
CP0=0.0
P0=10100
IF(IFLAG1)THEN
UPR=U0
ELSE
UPR=0.0
ENDIF
DO I=1,NXI
DO J=1,NYI
PR(I,J)=P0-0.5*(UTOT(I,J)**2-UPR**2)
ENDDO
ENDDO
RETURN
END