PROGRAM DRIVER REAL PI,GA,CL,CD,A,B COMMON/GRID/X(81,31),Y(81,31) COMMON/FLOW/PHI(81,31),RES(81,31),DELPHI(81,31) COMMON/PARAM/IMAX,JMAX,AMINF,DELT,COSANG,SINANG COMMON/STATS/NUMSUP,ITER,MAXRES,RMXI,RMXJ,MAXMACH, 1MMXI,MMXJ REAL RHOA1(81,31),RHOA3(81,31),RHOOUT(91) INTEGER RDES,RMXI,RMXJ,MMXI,MMXJ C-------------------------------------------------------------- C THIS PROGRAM SOLVES THE STEADY TRANSONIC FULL C EQUATION FOR FLOW PAST 2-D AIRFOILS C C-------------------------------------------------------------- Real GAMMA(31),CP(91) real MAXDEL,MAXRES,MAXMACH integer OCOUNT,ONUM PI = 4.*ATAN(1.) GA = 1.4 C read in grid data open(9,file='GRID.DAT',FORM='FORMATTED',STATUS='OLD') read(9,*) IMAX,JMAX read(9,400) ((Iloc,Jloc,X(I,J),Y(I,J),I=IMAX,1,-1),J=1,JMAX) WRITE(6,*) IMAX,JMAX close(9) 400 FORMAT(2I5,2F12.6) OPEN(UNIT=8,FILE='HISTORY') C C Get data from user C WRITE(6,*) 'Enter Mach number, ALPHA' READ(5,*) AMINF,ALPHA ALPHA=ALPHA*PI/180. COSANG=COS(ALPHA) SINANG=SIN(ALPHA) WRITE (6,*) 'AMINF=', AMINF , 1 'ALPHA=' , ALPHA * 180./PI WRITE (6,*) 'COSANG=' , COSANG WRITE (6,*) 'SINANG=' , SINANG WRITE(6,*) 'Restart? (1=y)' READ(5,*) RDES C C--Get data from user regarding to flap oscillation C C--------------------------------------------------------------------- C C Initialize PHI to zero or get restart data C C-------------------------------------------------------------------- IF (RDES .NE. 1) THEN DO J=1,JMAX GAMMA(J)= 0.0 DO I=1,IMAX PHI(I,J)=0.0 DELPHI(I,J) = 0.0 END DO END DO ELSE OPEN(UNIT=2,FILE='FPE_RESTART') WRITE(6,*) 'Reading restart file....' READ(2,*) AMINF,COSANG,SINANG WRITE(6,*) 'Mach number= ',AMINF WRITE(6,*) 'ALPHA= ',ACOS(COSANG)*180./PI READ(2,*) ((PHI(I,J),I=1,IMAX),J=1,JMAX) CLOSE(2) END IF 42 WRITE(6,*) 'Number of iterations, output increment?' READ(5,*) ITMAX , ONUM C C Main Iteration begins here C LCOUNT = 0 DO 100 ITER=1,ITMAX LCOUNT = LCOUNT + 1 IF (LCOUNT.GT.5) LCOUNT = 1 MAXDEL=0. MAXMACH=0. MAXRES=0. DELT = 0.1* (0.5 / 0.1) ** (0.25* FLOAT(LCOUNT-1)) CALL RESI(RHOA1,RHOA3,RHOOUT) CALL SOLVE(RHOA1,RHOA3,MAXDEL) C C Compute Gamma at TE C write (6,*) phi(imax-2,2), phi(3,2) GAMMA(1) = PHI(IMAX-2,2) - PHI(3,2) GAMMA(2) = GAMMA(1) WRITE (6,*)'iter numsup maxres rmxi rmxj MAXDEL gamma' WRITE (6,'(2I5,E14.8,2I5,E14.8,F10.4)') 1ITER,NUMSUP,MAXRES, 1RMXI,RMXJ,MAXDEL,GAMMA(1) * Print out iteration data IF(ITER/ONUM*ONUM.EQ.ITER) THEN CALL LOADS(RHOOUT,AMINF,GA,IMAX,X,Y,COSANG,SINANG, 1CP,CL,CD,CMC4,CMHNGE) c WRITE (6,'(I5,2F10.4)') (I,X(I,2),CP(I),I=2,IMAX-1) WRITE(6,*) 'CL=',CL WRITE(6,*) 'CD=',CD WRITE(6,*) 'CMC4=',CMC4 WRITE(8,*) ITER , CL , CD , CMC4 , CMHNGE END IF C Update top boundary DO 70 I=2,(IMAX/2) A=-ATAN2(((1.-AMINF**2.)**.5)*Y(I,JMAX),X(I,JMAX)) PHI(I,JMAX)=(Gamma(1)/(2.*PI))*A 70 CONTINUE DO 71 I=IMAX/2+1,IMAX-1 B=2*PI-ATAN2(((1.-AMINF**2.)**.5)*Y(I,JMAX),X(I,JMAX)) c WRITE(*,*) i,B PHI(I,JMAX)=(Gamma(1)/(2.*PI))*B 71 CONTINUE * Update left and right boundaries DO 80 J=2,JMAX PHI(1,J)=PHI(IMAX-2,J)-Gamma(1) PHI(IMAX,J)=PHI(3,J)+Gamma(1) 80 CONTINUE * Do bottom boundary DO 90 I=2,IMAX-1 C Metrics at I,J: xxs=(X(I+1,2)-X(I-1,2))/2. yxs=(Y(I+1,2)-Y(I-1,2))/2. xet=(X(I,3)-X(I,1))/2. yet=(Y(I,3)-Y(I,1))/2. yacL=1./(xxs*yet-xet*yxs) xsx=yet*yacL xsy=-xet*yacL etx=-yxs*yacL ety=xxs*yacL C--------------------------------------------------------------------- C C COMPUTE PHI(I,1): C C-------------------------------------------------------------------- PHIxs=(PHI(I+1,2)-PHI(I-1,2))/2. At=xsx*etx+xsy*ety Bt=etx*COSANG+ety*(SINANG) Ct=etx**2.+ety**2. c-------Ft=Fcorr*yacL------------------------------------------------- PHIEt=-(PHIxs*At+Bt)/Ct PHI(I,1)=PHI(I,3)-2.*PHIEt 90 CONTINUE C-------------------------------------------------------------------- C Final update of ghost boundaries: PHI(1,1)=PHI(IMAX-2,1)-Gamma(1) PHI(IMAX,1)=PHI(3,1)+Gamma(1) TIME = TIME + DELT 100 CONTINUE WRITE(*,*) 'Writing restart file' OPEN(UNIT=2,FILE='FPE_RESTART') WRITE(2,*) AMINF,COSANG,SINANG WRITE(2,*) ((PHI(I,J), I = 1 , IMAX) , J = 1 , JMAX) 110 CONTINUE WRITE(*,*) 'Writing CP file' OPEN(UNIT=3,FILE='CP_DATA') WRITE(3,121) '*' WRITE(3,*) 'Xsurf',CHAR(9),'Ysurf',CHAR(9),'CPsurf' DO 120 I=2,IMAX-1 WRITE(3,*) X(I,2),CHAR(9),Y(I,2),CHAR(9),(CP(I)+CP(I-1))/2 120 CONTINUE 121 FORMAT(A1) CLOSE(3) CLOSE(2) WRITE(6,*) 'Continue iteration? (1=y)' READ(5,*) I IF (I .EQ. 1) GO TO 42 STOP END SUBROUTINE SOLVE(RHOA1,RHOA3,MAXDEL) COMMON/FLOW/PHI(81,31),RES(81,31),DELPHI(81,31) COMMON/PARAM/IMAX,JMAX,AMINF,DELT,COSANG,SINANG INTEGER IMAX,JMAX REAL RHOA1(81,31),RHOA3(81,31) REAL A(91), B(91),C(91),RHS(91),TEMP(81,31) REAL MAXDEL C C AF2 Scheme C MAXDEL = 0.0 GA = 1.4 OMEGA = 1.5 C C XI - SWEEP C DO 1000 J = JMAX - 1 , 2 , -1 C C ASSEMBLE TRI-DIAGONAL MATRIX COEFFICIENTS C DO 10 I = 2 , IMAX - 1 B(I) = -1./DELT - (RHOA1(I,J)+RHOA1(I-1,J)) A(I) = RHOA1(I-1,J) C(I) = RHOA1(I,J) RHS(I) = OMEGA * RES(I,J)- 1./DELT * DELPHI(I,J+1) 10 CONTINUE C C PERFORM THOMAS ALGORITHM CALL THOMAS(A,B,C,2,IMAX-1,RHS) DO 20 I = 2 , IMAX - 1 20 DELPHI(I,J) = RHS(I) 1000 CONTINUE C C ETA SWEEP C DO 2000 I = 2 , IMAX - 1 DO 30 J = 2 , JMAX - 1 B(J) = -1./DELT -(RHOA3(I,J)+RHOA3(I,J-1)) 1 A(J) = RHOA3(I,j-1) C(J) = RHOA3(I,J) +1./DELT RHS(J) = -(DELPHI(I,J+1) - DELPHI(I,J))/DELT 30 CONTINUE CALL THOMAS(A,B,C,2,JMAX-1,RHS) DO 2000 J=2,JMAX-1 DELPHI(I,J) = RHS(J) PHI(I,J) = PHI(I,J) + DELPHI(I,J) MAXDEL = AMAX1(MAXDEL,ABS(DELPHI(I,J)/DELT)) 2000 CONTINUE C write (6,*) 'MAXDEL=', MAXDEL RETURN END SUBROUTINE THOMAS(A,B,C,ISTART,IEND,RHS) REAL A(91),B(91),C(91),X(91),RHS(91) C THOMAS ALGORITHM. IT IS ASSUMED ISTART < IEND C C RHS(ISTART-1) = 0.0 RHS(IEND+1) = 0.0 X(ISTART-1) = 0.0 A(ISTART) = 0.0 C(IEND) = 0.0 DO 10 I = ISTART , IEND DD = B(I) - A(I) * X(I-1) X(I) = C(I) / DD RHS(I) = (RHS(I) - A(I) * RHS(I-1)) / DD 10 CONTINUE DO 20 I = IEND , ISTART , -1 20 RHS(I) = RHS(I) - X(I) * RHS(I+1) RETURN END SUBROUTINE LOADS(RHOOUT,AMINF,GA,IMAX,X,Y,COSANG,SINANG, *CP,CL,CD,CMC4,CMHNGE) C REAL CX,CY,CL,CD,COSANG,SINANG,GA,AMINF REAL RHOOUT(91),X(81,31),Y(81,31),CP(91) C CL = 0.0 CD = 0.0 CMC4= 0.0 CMHNGE= 0.0 C DO 100 I = 1,IMAX CP(I) = 2.*(RHOOUT(I)**GA-1.)/(GA*AMINF**2.) CX = CP(I)*(X(I,2)-X(I+1,2)) CY = CP(I)*(Y(I+1,2)-Y(I,2)) CL = CL + CX*COSANG - CY*SINANG CD = CD + CX*SINANG + CY*COSANG XAVG = 0.5*(X(I,2)+X(I+1,2)) YAVG = 0.5*(Y(I,2)+Y(I+1,2)) CMC4 = CMC4 - CX*(XAVG-0.25) + CY*YAVG 100 CONTINUE C RETURN END SUBROUTINE RESI(RHOA1,RHOA3,RHOOUT) COMMON/PARAM/IMAX,JMAX,AMINF,DELT,COSANG,SINANG REAL RHOOUT(91) COMMON/GRID/X(81,31),Y(81,31) COMMON/FLOW/PHI(81,31),RES(81,31),DELPHI(81,31) COMMON/STATS/NUMSUP,ITER,MAXRES,RMXI,RMXJ,MAXMACH, 1MMXI,MMXJ INTEGER RDES,RMXI,RMXJ,MMXI,MMXJ C C Dimensions & counters: C C Quantities passed into the subroutine: C Quantities used in calculations: REAL RHO(100), $Ucon(100),Vcont(100),MsqLoc(100),Mu,RHOBar(100),RHOHat(100), $RHOUqJ(100),MAXRES,A,B,C,MAXMACH C C Quantities passed out of the subroutine: C REAL RHOA1(81,31),RHOA3(81,31) NUMSUP=0 GA = 1.4 DO 60 J=2,JMAX-1 C C Compute metrics, flow properties and UCON at (I+1/2,J) C DO 20 I=2,IMAX-1 XX = X(I+1,J) - X(I,J) XY = 0.25 * (X(I+1,J+1)-X(I+1,J-1)+X(I,J+1)-X(I,J-1)) YX = Y(I+1,J) - Y(I,J) YY = 0.25 * (Y(I+1,J+1)-Y(I+1,J-1)+Y(I,J+1)-Y(I,J-1)) YAC1 = 1. / (XX * YY - XY * YX) AX = YY * YAC1 AY = - XY * YAC1 BX = - YX * YAC1 BY = XX * YAC1 PX = PHI(I+1,J)-PHI(I,J) PY =.25*(PHI(I+1,J+1)-PHI(I+1,J-1)+PHI(I,J+1)-PHI(I,J-1)) IF (J .EQ. 2) THEN PY = -(PX*(AX *BX + AY * BY)+ COSANG * BX + SINANG * BY)/ $ (BX * BX + BY * BY) END IF c Calculate the u's and v's U = PX * AX + PY * BX+COSANG V = PX * AY + PY * BY+SINANG C Calculate density B=(1./(GA-1.)) A=1. - U * U -V * V C=(1+((GA-1.)/2)*(AMINF**2.)*A) RHO(I)=C**B C Calculate contravariant velocity over Jacobian Ucon(I) = (U * AX + V * AY) / YAC1 RHOA1(I,J) = (AX * AX + AY * AY) / YAC1 C Calculate local Mach number squared MsqLoc(I)=(u*u+v*v)/C IF (MsqLoc(I) .LT. 0.) THEN WRITE(*,*) 'uh oh... Msq<0' STOP ELSE IF (SQRT(MsqLoc(I)) .GT. MAXMACH) THEN MAXMACH=SQRT(MsqLoc(I)) MMXI=I MMXJ=J END IF END IF C Supersonic counter IF (MsqLoc(I) .GT. 1.) THEN NumSup=NumSup+1 END IF 20 CONTINUE C Set ghost boundaries MsqLoc(1)=MsqLoc(IMAX-2) MsqLoc(IMAX)=MsqLoc(3) RHO(1)=RHO(IMAX-2) RHO(IMAX)=RHO(3) c WRITE(*,*) 'RHO(',IMAX-1,')=',RHO(IMAX-1), c $'RHO(2)= ',RHO(2) C Loop to calculate RHObar DO 30 I=2,IMAX-1 IF (Ucon(I) .GE. 0.) THEN Mu=MAX(0.,1.-2./(MsqLoc(I-1)+MsqLoc(I))) RHOBar(I)=RHO(I)-Mu*(RHO(I)-RHO(I-1)) ELSE Mu=MAX(0.,1.-2./(MsqLoc(I)+MsqLoc(I+1))) RHOBar(I)=RHO(I)-Mu*(RHO(I)-RHO(I+1)) END IF 30 CONTINUE C Calculate quantities (RHOBar*U/J) and (RHOBar*A1) DO 40 I=2,IMAX-1 RHOUqJ(I)=RHOBar(I)*Ucon(I) RHOA1(I,J)=RHOBar(I) * RHOA1(I,J) 40 CONTINUE c Set ghost boundaries RHOUqJ(1)=RHOUqJ(IMAX-2) RHOA1(1,J)=RHOA1(IMAX-2,J) DO 50 I=2,IMAX-1 Res(I,J)=RHOUqJ(I)-RHOUqJ(I-1) 50 CONTINUE IF (J .EQ. 2) THEN DO 55 I=2,IMAX-1 RHOOUT(I) = RHO(I) 55 CONTINUE RHOOUT(IMAX) = RHO(3) RHOOUT(1) = RHO(IMAX-2) END IF 60 CONTINUE C Now repeat process going from bottom to top, calculating C the J+1/2 half-point quantities, then move over a line and repeat. DO 100 I=2,IMAX-1 DO 70 J=2,JMAX-1 C Calculation of the PHI-derivatives wrt Xsi and ETA: PX = .25*(PHI(I+1,J+1)-PHI(I-1,J+1)+PHI(I+1,J)-PHI(I-1,J)) PY = PHI(I,J+1) - PHI(I,J) XX = 0.25 * (X(I+1,J+1)-X(I-1,J+1)+X(I+1,J)-X(I-1,J)) YX = 0.25 * (Y(I+1,J+1)-Y(I-1,J+1)+Y(I+1,J)-Y(I-1,J)) XY = X(I,J+1) - X(I,J) YY = Y(I,J+1) - Y(I,J) YAC1 = 1. / (XX * YY - XY * YX) AX = YY * YAC1 AY = - XY * YAC1 BX = - YX * YAC1 BY = XX * YAC1 C Calculate the u's and v's u = PX * AX+PY * BX +COSANG v = PX * AY+PY * BY+SINANG C Calculate density, called RHOHat here B=(1./(GA-1.)) A=(1.-u**2.-v**2 ) C=(1+((GA-1.)/2)*(AMINF**2.)*A) RHOHat(J)=C**B c Calculate contravariant velocities Vcont(J)= RHOHAT(J) * ( U * BX + V * BY) / YAC1 RHOA3(I,J) = RHOHAT(J) * (BX * BX + BY * BY) / YAC1 70 CONTINUE c--------------------------------------------------------------------- C Add the contribution to Residual(I,J) C (Treat J=2 differently here) Res(I,2)=Res(I,2)+(2.*VCONT(2)) IF (ABS(Res(I,2)) .GT. MAXRES) THEN MAXRES=ABS(Res(I,2)) RMXI=I RMXJ=2 END IF C Loop from J=3 to J=JMAX-1 DO 80 J=3,JMAX-1 Res(I,J)=Res(I,J)+ ( VCONT(J) - VCONT(J-1)) IF (ABS(Res(I,J)) .GT. MAXRES) THEN MAXRES=ABS(Res(I,J)) RMXI=I RMXJ=J END IF 80 CONTINUE 100 CONTINUE write (6,*) 'MAXRES=', MAXRES RETURN END