c************************************************************************ c * c PROGRAM THINFOIL * c * c FINITE-DIFFERENCE SOLUTION OF INCOMPRESSIBLE POTENTIAL FLOW * c PAST A THIN SYMMETRIC AIRFOIL AT ZERO ANGLE OF ATTACK * c * c SOR, SLOR and ADI Iterative schemes * c * c by Valery Razgonyaev for W.H. Mason * c * c inspired by MORAN's THINAIR code, PAGE 326-330 * c * c code emphasizes clarity over efficiency (hopefully) * c * c Last modification: April 9, 1992 * c************************************************************************ IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*1 ICASE COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91), 1 X(91),Y(41),NX,NY,ILE,ITE, 2 AX(91),BX(91),CX(91),AY(41),BY(41),CY(41), 3 RELAX,TCMAX DATA IREAD /5/ IWRIT /6/ IFILE /7/ OPEN(UNIT = IFILE, FILE = 'THINAIR.OUT', STATUS = 'UNKNOWN') WRITE(IWRIT,1000) WRITE(IFILE,1000) c define the mesh, initialize the potential, define bc's call setup c solve the equations with iteration 10 call solve c output the solution on the surface call output c ask if more iterations are desired write(IWRIT,1040) write(IFILE,1040) read(IREAD,'(a)') icase write(IFILE,1050) icase KOUNT = KOUNT+maxit IF (icase .eq. 'Y' .or. ICASE .eq. 'y') go to 10 close(ifile) stop c formats 1000 format(/5X,'Program THINFOIL'//5x,'Finite difference solution of'/ 1 5x,'the Laplace Eqn. over a biconvex airfoil'/) 1040 format(/5x,'Additional iteration? (Y/N):',$) 1050 format(5x,A1/) END subroutine solve c select an iteration scheme and solve the equations IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*1 ICASE CHARACTER*4 SCHEME(4) COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91), 1 X(91),Y(41),NX,NY,ILE,ITE, 2 AX(91),BX(91),CX(91),AY(41),BY(41),CY(41), 3 RELAX,TCMAX DATA IREAD /5/ IWRIT /6/ IFILE /7/ ************************************************************************* * INPUT ITERATIVE SCHEME TYPE (ITYPE), * * MAXIMUM NUMBER OF ITERATIONS (maxit) * * AND OVERRELAXATION FACTOR (RELAX) * ************************************************************************* SCHEME(1) = 'SOR ' SCHEME(2) = 'SLOR' SCHEME(3) = 'AF-1' SCHEME(4) = 'AF-2' KOUNT = 0 030 WRITE(IWRIT,1000) READ(IREAD,*) ITYPE WRITE(IWRIT,1005) SCHEME(ITYPE) WRITE(IFILE,1005) SCHEME(ITYPE) WRITE(IWRIT,1010) READ(IREAD,*) RELAX write(iwrit,1015) read(iread,*) maxit WRITE(IFILE,1020) maxit, RELAX WRITE(IWRIT,1030) WRITE(IFILE,1030) ************************************************************************* * COMPUTATION OF THE POTENTIAL * ************************************************************************* DO 80 iter = 1,maxit ************************************************************************* * SET UPDATE OUTSIDE THE COMPUTATIONAL SPACE * ************************************************************************* DO 40 I = 2,NX C(I,1) = C(I,3) C(I,NY+1) = C(I,NY-1) 40 continue DO 50 J = 2,NY C(1,J) = C(3,J) C(NX+1,J) = C(NX-1,J) 50 continue ************************************************************************* * UPDATE THE POTENTIAL * ************************************************************************* DO 60 I = 1,NX+1 DO 60 J = 1,NY+1 PHI(I,J) = PHI(I,J)+C(I,J) 60 continue ************************************************************************* * UPDATE THE RESIDUAL * ************************************************************************* do 70 I = 2,NX do 70 J = 2,NY R(I,J) = AX(I)*PHI(I-1,J)+BX(I)*PHI(I,J)+CX(I)*PHI(I+1,J)+ 1 AY(J)*PHI(I,J-1)+BY(J)*PHI(I,J)+CY(J)*PHI(I,J+1) 70 continue ************************************************************************* * OUTPUT ITERATION HISTORY * ************************************************************************* CALL STAT(KOUNT+iter,itype) ************************************************************************* * COMPUTATION OF THE POTENTIAL UPDATE * ************************************************************************* IF (ITYPE .EQ. 1) THEN CALL SOR ELSE IF (ITYPE .EQ. 2) THEN CALL SLOR ELSE IF (ITYPE .EQ. 3) THEN CALL AF1(iter) ELSE CALL AF2(iter) END IF 80 continue return 1000 format(//5X,'Iterative Scheme Options:'// 1 8x,'1 - SOR '/ 2 8x,'2 - SLOR'/ 3 8x,'3 - AF-1'/ 4 8x,'4 - AF-2'/ 5 20x,'Select option:',$) 1005 format(/5X,'Iterative Scheme Selected: - ',A4/) 1010 format(/5X,'INPUT relaxation factor:'/ 1 8x,'Suggested values for the ADI Relaxation Factors:'/ 2' AF-1 - 2.0'/ 3' AF-2 - 1.333 :',$) 1015 format(//5X,'INPUT number of iterations, maxit:'$) 1020 format(/4X,I6,F20.3) 1030 format(//33x,'Convergence History'/ 1 20x,'Residual', 27x,'Correction'/ 2 9x,'-------------------------------', 3 6x,'-------------------------------'/ 4 3X,'ITER IMAX JMAX MAX AVRG', 5 7x, 'IMAX JMAX MAX AVRG') end subroutine setup c setup mesh, initialize th epotential and define bc's IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*1 ICASE CHARACTER*4 SCHEME(4) COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91), 1 X(91),Y(41),NX,NY,ILE,ITE, 2 AX(91),BX(91),CX(91),AY(41),BY(41),CY(41), 3 RELAX,TCMAX DATA IREAD /5/ IWRIT /6/ IFILE /7/ c SET UP MESH CALL MESH c MAKE THE INITIAL GUESS FOR POTENTIAL & UPDATE DO 10 I = 1,NX+1 DO 10 J = 1,NY+1 PHI(I,J) = 0.D0 C(I,J) = 0.D0 10 continue c DEFINE THE BOUNDARY CONDITION ON THE SURFACE * tcmax = 0.05D0 do 20 i = 1,nx fp(i) = 0.0 if(i .ge. ile .and. i .le. ite) then fp(i) = 2.D0*tcmax*(1.D0-2.D0*X(i)) endif PHI(i,1) =-(Y(3)-Y(1))*fp(i) 20 continue c output the grid and surface slope values WRITE(IFILE,4020) WRITE(Iwrit,4020) i = 1 write(iwrit,3031) i,x(i),fp(i) WRITE(IFILE,3030) I,X(I),fp(i) DO 80 I = 2,NX dxgrid = x(i) - x(i-1) write(iwrit,3030) i,x(i),dxgrid,fp(i) 80 WRITE(IFILE,3030) I,X(I),dxgrid,fp(i) WRITE(IFILE,4030) WRITE(iwrit,4030) j = 1 WRITE(Iwrit,3030) J,Y(J) WRITE(IFILE,3030) J,Y(J) DO 90 J = 2,NY dygrid = y(j) - y(j-1) WRITE(Iwrit,3030) J,Y(J),dygrid 90 WRITE(IFILE,3030) J,Y(J),dygrid return 3030 format(4x,I3,5E14.5) 3031 format(4x,i3,e14.5,14x,3e14.5) 4020 format(/8X,'X GRID and slope BC, fp'// 1 6x,'i',6x,'x(i)',7x,'x(i)-x(i-1)',6x,'fp(i)') 4030 format(/8X,'Y GRID'// 1 6x,'j',6x,'y(j)',7x,'y(j)-y(j-1)') 4040 format(/6X,I6,f20.3) end SUBROUTINE SOR ************************************************************** * POINT GAUSS-SEIDEL (SOR) ITERATIVE SCHEME * ************************************************************** IMPLICIT REAL*8 (A-H,O-Z) COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91), 1 X(91),Y(41),NX,NY,ILE,ITE, 2 AX(91),BX(91),CX(91),AY(41),BY(41),CY(41), 3 RELAX,TCMAX DO 100 I = 2,NX DO 100 J = 2,NY AXI = AX(I) AYJ = AY(J) IF (I. EQ. 2) 1 AXI = 0.D0 IF (I. EQ. NX) 1 AXI =-BX(NX) IF (J .EQ. 2) 1 AYJ = 0.D0 IF (J .EQ. NY) 1 AYJ =-BY(NY) C(I,J) =-RELAX*(R(I,J)+AXI*C(I-1,J)+AYJ*C(I,J-1))/ 1 (BX(I)+BY(J)) 100 continue RETURN END SUBROUTINE SLOR ************************************************************** * LINE GAUSS-SEIDEL (SLOR) ITERATIVE SCHEME * ************************************************************** IMPLICIT REAL*8 (A-H,O-Z) REAL*8 MAIND(41),LOWD(41),UPPD(41),RHS(41) COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91), 1 X(91),Y(41),NX,NY,ILE,ITE, 2 AX(91),BX(91),CX(91),AY(41),BY(41),CY(41), 3 RELAX,TCMAX ************************************************************************* * SWEEP ALONG Y DIRECTION * ************************************************************************* DO 210 I = 2,NX ************************************************************************* * SET TRIDIAGONAL MATRIX & RHS VECTOR FOR UPDATE IN I-th COLUMN * ************************************************************************* DO 200 J = 2,NY AXI = AX(I) AYJ = AY(J) CYJ = CY(J) IF (I. EQ. 2) 1 AXI = 0.D0 IF (I. EQ. NX) 1 AXI =-BX(NX) IF (J .EQ. 2) 1 CYJ =-BY(2) IF (J .EQ. NY) 1 AYJ =-BY(NY) MAIND(J-1) = BX(I)+BY(J) UPPD(J-1) = CYJ LOWD(J-1) = AYJ RHS(J-1) =-RELAX*(R(I,J)+AXI*C(I-1,J)) 200 continue ************************************************************************* * CALL FOR THOMAS ALGORITHM * ************************************************************************* CALL THOMAS(LOWD,MAIND,UPPD,RHS,NY-1) DO 210 J = 2,NY C(I,J) = RHS(J-1) 210 continue RETURN END SUBROUTINE AF1(NIT) ************************************************************************ * APPROXIMATE FACTORISATION AF-1 (ADI) ITERATIVE SCHEME * ************************************************************************ IMPLICIT REAL*8 (A-H,O-Z) REAL*8 F(91,41),MAIND(91),LOWD(91),UPPD(91),RHS(91),M COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91), 1 X(91),Y(41),NX,NY,ILE,ITE, 2 AX(91),BX(91),CX(91),AY(41),BY(41),CY(41), 3 RELAX,TCMAX K = NIT-1 M = 11.D0 P = 2.D0*((K-M*INT(K/M))/(M-1.D0)-1.D0) ALFA = (0.5D0*(Y(3)-Y(2)))**P ************************************************************************* * INVERT OUTER OPERATOR * * OR * * SWEEP ALONG I DIRECTION * ************************************************************************* DO 320 J = 2,NY DO 310 I = 2,NX AXI = AX(I) CXI = CX(I) IF (I .EQ. 2) 1 CXI =-BX(2) IF (I .EQ. NX) 1 AXI =-BX(NX) MAIND(I-1) = BX(I)-ALFA UPPD(I-1) = CXI LOWD(I-1) = AXI RHS(I-1) = ALFA*RELAX*R(I,J) 310 continue c************************************************************************ c CALL FOR THOMAS ALGORITHM * c************************************************************************ CALL THOMAS(LOWD,MAIND,UPPD,RHS,NX-1) DO 320 I = 2,NX F(I,J) = RHS(I-1) 320 continue c************************************************************************ c INVERT INNER OPERATOR * c OR * c SWEEP ALONG J DIRECTION * c************************************************************************ DO 340 I = 2,NX DO 330 J = 2,NY AYJ = AY(J) CYJ = CY(J) IF (J .EQ. 2) 1 CYJ =-BY(2) IF (J .EQ. NY) 1 AYJ =-BY(NY) MAIND(J-1) = BY(J)-ALFA UPPD(J-1) = CYJ LOWD(J-1) = AYJ RHS(J-1) = F(I,J) 330 continue c************************************************************************ c CALL FOR THOMAS ALGORITHM * c************************************************************************ CALL THOMAS(LOWD,MAIND,UPPD,RHS,NY-1) DO 340 J = 2,NY C(I,J) = RHS(J-1) 340 continue RETURN END SUBROUTINE AF2(NIT) c************************************************************************ c APPROXIMATE FACTORISATION AF-2 (ADI) ITERATIVE SCHEME * c************************************************************************ IMPLICIT REAL*8 (A-H,O-Z) REAL*8 F(91,41),MAIND(91),LOWD(91),UPPD(91),RHS(91),M COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91), 1 X(91),Y(41),NX,NY,ILE,ITE, 2 AX(91),BX(91),CX(91),AY(41),BY(41),CY(41), 3 RELAX,TCMAX K = NIT-1 M = 11.D0 P = (K-M*INT(K/M))/(M-1.D0)-1.D0 ALFA = ((Y(3)-Y(2))/2.D0)**P ************************************************************************* * INVERT OUTER OPERATOR * * OR * * SWEEP ALONG I DIRECTION * ************************************************************************* DO 410 J = 2,NY DO 400 I = 2,NX DLOW = 0.D0 IF (I .EQ. NX) 1 DLOW =-1.D0 MAIND(I-1) = ALFA*(X(I+1)-X(I))+1.D0 UPPD(I-1) =-1.D0 LOWD(I-1) = DLOW RHS(I-1) =-ALFA*RELAX*R(I,J)*(X(I+1)-X(I)) 400 continue ************************************************************************* * CALL FOR THOMAS ALGORITHM * ************************************************************************* CALL THOMAS(LOWD,MAIND,UPPD,RHS,NX-1) DO 410 I = 2,NX F(I,J) = RHS(I-1) 410 continue ************************************************************************* * INVERT INNER OPERATOR * * OR * * SWEEP ALONG J DIRECTION * ************************************************************************* DO 430 I = 2,NX DO 420 J = 2,NY AYJ = AY(J) CYJ = CY(J) CIJ = C(I-1,J) IF (J .EQ. 2) 1 CYJ =-BY(2) IF (J .EQ. NY) 1 AYJ =-BY(NY) IF (I .EQ. 2) 1 CIJ = 0.D0 MAIND(J-1) = BY(J)-ALFA/(X(I)-X(I-1)) UPPD(J-1) = CYJ LOWD(J-1) = AYJ RHS(J-1) = F(I,J)-ALFA*CIJ/(X(I)-X(I-1)) 420 continue ************************************************************************* * CALL FOR THOMAS ALGORITHM * ************************************************************************* CALL THOMAS(LOWD,MAIND,UPPD,RHS,NY-1) DO 430 J = 2,NY C(I,J) = RHS(J-1) 430 continue RETURN END SUBROUTINE THOMAS(A,B,C,D,N) c************************************************************************ c INVERT TRIDIAGONAL MATRIX USING THOMAS ALGORITHM * c************************************************************************ REAL*8 A(91), B(91), C(91), D(91), M DO 400 K = 2,N M = A(K)/B(K-1) B(K) = B(K)-M*C(K-1) D(K) = D(K)-M*D(K-1) 400 continue D(N) = D(N)/B(N) DO 410 K = N-1,1,-1 D(K) = (D(K)-C(K)*D(K+1))/B(K) 410 continue RETURN END SUBROUTINE MESH c*************************************************** c MESH GENERATION * c*************************************************** IMPLICIT REAL*8 (A-H,O-Z) COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91), 1 X(91),Y(41),NX,NY,ILE,ITE, 2 AX(91),BX(91),CX(91),AY(41),BY(41),CY(41), 3 RELAX,TCMAX DATA PI /3.141592685/ IREAD /5/ IWRIT /6/ IFILE /7/ c************************************************************************ c INPUT FIELD DIMENSIONS (XW < X < XE, 0 < Y < YN) * c AND NUMBERS OF INTERVALS UPSTREAM (NUP) * c DOWNSTREAM (NDOWN), ON(NON), AND ABOVE (NABOVE) * c THE AIRFOIL * c************************************************************************ WRITE(IWRIT,2010) READ(IREAD,*) XW,XE,YN WRITE(IWRIT,2020) READ(IREAD,*) NUP WRITE(IWRIT,2022) READ(IREAD,*) NDOWN WRITE(IWRIT,2024) READ(IREAD,*) NON WRITE(IWRIT,2026) READ(IREAD,*) NABOVE ILE = NUP + 2 ITE = ILE + NON NX = ITE + NDOWN NY = NABOVE + 2 IF (NX .GT. 86 ) THEN WRITE(IWRIT,2030) NX STOP END IF IF (NY .GT. 36) THEN WRITE(IWRIT,2040) NY STOP END IF WRITE(IFILE,2010) WRITE(IFILE,2050) XW,XE,YN WRITE(IFILE,2020) WRITE(IFILE,2060) NUP,NDOWN,NON,NABOVE c************************************************************************ c GENERATE MESH * c************************************************************************ DO 510 I = 2,ILE 510 X(I) = XW*(1.0D0-DCOS(PI/2.D0*(ILE-I)/DBLE(ILE-2))) X(1) = 2.D0*X(2)-X(3) DO 520 I = ILE,ITE 520 X(I) = 0.5D0*(1.D0-COS(PI*(I-ILE)/DBLE(ITE-ILE))) DO 530 I = ITE,NX 530 X(I) = XE-(XE-1.D0)*COS(PI/2.D0*(I-ITE)/DBLE(NX-ITE)) X(NX+1) = X(NX)+(X(NX)-X(NX-1)) DO 540 J = 2,NY 540 Y(J) = YN*(1.D0-COS(PI/2.D0*(J-2)/DBLE(NY-2))) Y(1) = 2.D0*Y(2) - Y(3) Y(NY+1) = Y(NY) + (Y(NY) - Y(NY-1)) c************************************************************************ c COMPUTE COEFFICIENTS FOR LAPLACIAN * c * c L(PHI(I,J)) = AX(I)*PHI(I-1,J)+BX(I)*PHI(I,J)+CX(I)*PHI(I+1,J)+ * c * c AY(I)*PHI(I,J-1)+BY(J)*PHI(I,J)+CY(J)*PHI(I,J+1) * c * c************************************************************************ DO 550 I = 2,NX AX(I) = 2.D0/(X(I+1)-X(I-1))/(X(I)-X(I-1)) CX(I) = 2.D0/(X(I+1)-X(I-1))/(X(I+1)-X(I)) BX(I) =-(AX(I)+CX(I)) 550 continue DO 560 J = 2,NY AY(J) = 2.D0/(Y(J+1)-Y(J-1))/(Y(J)-Y(J-1)) CY(J) = 2.D0/(Y(J+1)-Y(J-1))/(Y(J+1)-Y(J)) BY(J) =-(AY(J)+CY(J)) 560 continue RETURN 2010 format(/5X, 1 'INPUT grid boundary locations; Xmin, Xmax, Ymax:',$) 2020 format(/5X, 1 'INPUT number of points upstream of airfoil, Nup:',$) 2022 format(5x, 1 ' downstream of airfoil, Ndown:',$) 2024 format(5X, 1 ' on the airfoil, NON:',$) 2026 format(5X, 1 ' above the airfoil, NABOVE:',$) 2030 format(/5x,'Too many x grid points. Max is 86, your nx is',i3/) 2040 format(/5x,'Too many y grid points. Max is 36, your ny is',i3/) 2050 format(/4X,3F12.4) 2060 format(/4X,4I6) END SUBROUTINE STAT(KOUNT,itype) c***************************************************** c Solution HISTORY * c***************************************************** IMPLICIT REAL*8 (A-H,O-Z) COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91), 1 X(91),Y(41),NX,NY,ILE,ITE, 2 AX(91),BX(91),CX(91),AY(41),BY(41),CY(41), 3 RELAX,TCMAX DATA IREAD /5/ IWRIT /6/ IFILE /7/ c************************************************************************ c FIND MAXIMAL RESIDUAL AND UPDATE * c************************************************************************ IMAXR = 2 JMAXR = 2 IMAXC = 2 JMAXC = 2 RAVRG = 0.D0 CAVRG = 0.D0 DO 600 I = 2,NX DO 600 J = 2,NY IF (DABS(R(I,J)) .GT. DABS(R(IMAXR,JMAXR))) THEN IMAXR = I JMAXR = J END IF IF (DABS(C(I,J)) .GT. DABS(C(IMAXC,JMAXC))) THEN IMAXC = I JMAXC = J END IF RAVRG = RAVRG+R(I,J)**2 CAVRG = CAVRG+C(I,J)**2 600 continue RMAX = DABS(R(IMAXR,JMAXR)) CMAX = DABS(C(IMAXC,JMAXC)) RAVRG = DSQRT(RAVRG/DBLE((NX-1)*(NY-1))) CAVRG = DSQRT(CAVRG/DBLE((NX-1)*(NY-1))) c if(itype .lt. 3 .and. kount*(kount/10) .ne. kount) return c************************************************************************ c OUTPUT ITERATION HISTORY * c************************************************************************ WRITE(IFILE,3010) KOUNT,IMAXR,JMAXR,RMAX, 1 RAVRG,IMAXC,JMAXC,CMAX,CAVRG WRITE(IWRIT,3010) KOUNT,IMAXR,JMAXR,RMAX, 1 RAVRG,IMAXC,JMAXC,CMAX,CAVRG RETURN 3000 format(/3X,A50) 3010 format(1X,3I5,2E12.4,3X,2I5,2E12.4) END SUBROUTINE output c************************************************************************ c COMPUTE AND PRINT PRESSURE DISTRIBUTION ON THE AIRFOIL * c************************************************************************ IMPLICIT REAL*8 (A-H,O-Z) COMMON /SOL/ PHI(91,41),C(91,41),R(91,41),FP(91), 1 X(91),Y(41),NX,NY,ILE,ITE, 2 AX(91),BX(91),CX(91),AY(41),BY(41),CY(41), 3 RELAX,TCMAX DATA IREAD /5/ IWRIT /6/ IFILE /7/ WRITE(IWRIT,4000) WRITE(IFILE,4000) ILEP = ILE + 1 DO 700 I = ILEP,ITE XP = 0.5D0*(X(I) + X(I-1)) U = (PHI(I,2) - PHI(I-1,2))/(X(I) - X(I-1)) CP = - 2.D0*U WRITE(IWRIT,4010) I,XP,U,CP WRITE(IFILE,4010) I,XP,U,CP 700 continue 4000 format(//3X,' PRESSURE DISTRIBUTION ON THE AIRFOIL'//5X, 1 'I',10X,'x/c',13X,'U',14X,'Cp') 4010 format(3X,I3,3F15.5) RETURN END