C*********************************************************************** C Module: xpanel.f C C Copyright (C) 2000 Mark Drela C C This program is free software; you can redistribute it and/or modify C it under the terms of the GNU General Public License as published by C the Free Software Foundation; either version 2 of the License, or C (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public License C along with this program; if not, write to the Free Software C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. C*********************************************************************** SUBROUTINE APCALC INCLUDE 'XFOIL.INC' C C---- set angles of airfoil panels DO 10 I=1, N-1 SX = X(I+1) - X(I) SY = Y(I+1) - Y(I) IF(SX.EQ.0.0 .AND. SY.EQ.0.0) THEN IF(I.EQ.1) THEN APANEL(I) = ATAN2( -NY(I) , -NX(I) ) ELSE APANEL(I) = ATANC( -NY(I) , -NX(I) , APANEL(I-1) ) ENDIF ELSE IF(I.EQ.1) THEN APANEL(I) = ATAN2( SX , -SY ) ELSE APANEL(I) = ATANC( SX , -SY , APANEL(I-1) ) ENDIF ENDIF 10 CONTINUE C C---- TE panel I = N IP = 1 IF(SHARP) THEN APANEL(I) = PI ELSE SX = X(IP) - X(I) SY = Y(IP) - Y(I) ccc APANEL(I) = ATAN2( -SX , SY ) + PI APANEL(I) = ATANC( SX , -SY , APANEL(I-1) ) ENDIF C RETURN END SUBROUTINE NCALC(X,Y,S,N,XN,YN) C--------------------------------------- C Calculates normal unit vector C components at airfoil panel nodes C--------------------------------------- DIMENSION X(N), Y(N), S(N), XN(N), YN(N) C IF(N.LE.1) RETURN C CALL SEGSPL(X,XN,S,N) CALL SEGSPL(Y,YN,S,N) DO 10 I=1, N SX = YN(I) SY = -XN(I) SMOD = SQRT(SX*SX + SY*SY) XN(I) = SX/SMOD YN(I) = SY/SMOD 10 CONTINUE C C---- average normal vectors at corner points DO 20 I=1, N-1 IF(S(I) .EQ. S(I+1)) THEN SX = 0.5*(XN(I) + XN(I+1)) SY = 0.5*(YN(I) + YN(I+1)) SMOD = SQRT(SX*SX + SY*SY) XN(I) = SX/SMOD YN(I) = SY/SMOD XN(I+1) = SX/SMOD YN(I+1) = SY/SMOD ENDIF 20 CONTINUE C RETURN END SUBROUTINE PSILIN(I,XI,YI,NXI,NYI,PSI,PSI_NI,GEOLIN,SIGLIN) C----------------------------------------------------------------------- C Calculates current streamfunction Psi at panel node or wake node C I due to freestream and all bound vorticity Gam on the airfoil. C Sensitivities of Psi with respect to alpha (Z_ALFA) and inverse C Qspec DOFs (Z_QDOF0,Z_QDOF1) which influence Gam in inverse cases. C Also calculates the sensitivity vector dPsi/dGam (DZDG). C C If SIGLIN=True, then Psi includes the effects of the viscous C source distribution Sig and the sensitivity vector dPsi/dSig C (DZDM) is calculated. C C If GEOLIN=True, then the geometric sensitivity vector dPsi/dn C is calculated, where n is the normal motion of the jth node. C C Airfoil: 1 < I < N C Wake: N+1 < I < N+NW C----------------------------------------------------------------------- INCLUDE 'XFOIL.INC' REAL NXO, NYO, NXP, NYP, NXI, NYI LOGICAL GEOLIN,SIGLIN REAL THETJ(0:IZX) C C---- distance tolerance for determining if two points are the same SEPS = (S(N)-S(1)) * 1.0E-5 C IO = I C COSA = COS(ALFA) SINA = SIN(ALFA) C DO JO=1, N DZDG(JO) = 0.0 DZDN(JO) = 0.0 DQDG(JO) = 0.0 ENDDO C DO JO=1, N DZDM(JO) = 0.0 DQDM(JO) = 0.0 ENDDO C Z_QINF = 0. Z_ALFA = 0. Z_QDOF0 = 0. Z_QDOF1 = 0. Z_QDOF2 = 0. Z_QDOF3 = 0. C PSI = 0. PSI_NI = 0. C QTAN1 = 0. QTAN2 = 0. QTANM = 0. C IF(SHARP) THEN SCS = 1.0 SDS = 0.0 ELSE SCS = ANTE/DSTE SDS = ASTE/DSTE ENDIF C THETJ(0) = 0.0 C DO 10 JO=1, N JP = JO+1 C JM = JO-1 JQ = JP+1 C IF(JO.EQ.1) THEN JM = JO ELSE IF(JO.EQ.N-1) THEN JQ = JP ELSE IF(JO.EQ.N) THEN JP = 1 IF((X(JO)-X(JP))**2 + (Y(JO)-Y(JP))**2 .LT. SEPS**2) GO TO 12 ENDIF C DSO = SQRT((X(JO)-X(JP))**2 + (Y(JO)-Y(JP))**2) C C------ skip null panel IF(DSO .EQ. 0.0) GO TO 10 C DSIO = 1.0 / DSO C APAN = APANEL(JO) C RX1 = XI - X(JO) RY1 = YI - Y(JO) RX2 = XI - X(JP) RY2 = YI - Y(JP) C SX = (X(JP) - X(JO)) * DSIO SY = (Y(JP) - Y(JO)) * DSIO C X1 = SX*RX1 + SY*RY1 X2 = SX*RX2 + SY*RY2 YY = SX*RY1 - SY*RX1 C RS1 = RX1*RX1 + RY1*RY1 RS2 = RX2*RX2 + RY2*RY2 C C------ set log(r^2) and arctan(x/y), correcting for reflection if any IF(IO.NE.JO .AND. RS1.GT.0.0) THEN G1 = LOG(RS1) ccc T1 = ATAN2(X1,YY) T1 = ATANC(X1,YY,THETJ(JO-1)) ELSE G1 = 0.0 T1 = 0.0 ENDIF THETJ(JO) = T1 C IF(IO.NE.JP .AND. RS2.GT.0.0) THEN G2 = LOG(RS2) ccc T2 = ATAN2(X2,YY) T2 = ATANC(X2,YY,THETJ(JO)) ELSE G2 = 0.0 T2 = 0.0 ENDIF THETJ(JP) = T2 C X1I = SX*NXI + SY*NYI X2I = SX*NXI + SY*NYI YYI = SX*NYI - SY*NXI C IF(GEOLIN) THEN NXO = NX(JO) NYO = NY(JO) NXP = NX(JP) NYP = NY(JP) C X1O =-((RX1-X1*SX)*NXO + (RY1-X1*SY)*NYO)*DSIO-(SX*NXO+SY*NYO) X1P = ((RX1-X1*SX)*NXP + (RY1-X1*SY)*NYP)*DSIO X2O =-((RX2-X2*SX)*NXO + (RY2-X2*SY)*NYO)*DSIO X2P = ((RX2-X2*SX)*NXP + (RY2-X2*SY)*NYP)*DSIO-(SX*NXP+SY*NYP) YYO = ((RX1+X1*SY)*NYO - (RY1-X1*SX)*NXO)*DSIO-(SX*NYO-SY*NXO) YYP =-((RX1-X1*SY)*NYP - (RY1+X1*SX)*NXP)*DSIO ENDIF C IF(JO.EQ.N) GO TO 11 C IF(SIGLIN) THEN C C------- set up midpoint quantities X0 = 0.5*(X1+X2) RS0 = X0*X0 + YY*YY G0 = LOG(RS0) ccc T0 = ATAN2(X0,YY) T0 = ATANC(X0,YY,THETJ(JO)) C C------- calculate source contribution to Psi for 1-0 half-panel DXINV = 1.0/(X1-X0) PSUM = X0*(T0-APAN) - X1*(T1-APAN) + 0.5*YY*(G1-G0) PDIF = ((X1+X0)*PSUM + RS1*(T1-APAN) - RS0*(T0-APAN) & + (X0-X1)*YY) * DXINV C PSX1 = -(T1-APAN) PSX0 = T0-APAN PSYY = 0.5*(G1-G0) C PDX1 = ((X1+X0)*PSX1 + PSUM + 2.0*X1*(T1-APAN) - PDIF) * DXINV PDX0 = ((X1+X0)*PSX0 + PSUM - 2.0*X0*(T0-APAN) + PDIF) * DXINV PDYY = ((X1+X0)*PSYY + 2.0*(X0-X1 + YY*(T1-T0)) ) * DXINV C DSM = SQRT((X(JP)-X(JM))**2 + (Y(JP)-Y(JM))**2) DSIM = 1.0/DSM C CCC SIG0 = (SIG(JP) - SIG(JO))*DSIO CCC SIG1 = (SIG(JP) - SIG(JM))*DSIM CCC SSUM = SIG0 + SIG1 CCC SDIF = SIG0 - SIG1 C SSUM = (SIG(JP) - SIG(JO))*DSIO + (SIG(JP) - SIG(JM))*DSIM SDIF = (SIG(JP) - SIG(JO))*DSIO - (SIG(JP) - SIG(JM))*DSIM C PSI = PSI + QOPI*(PSUM*SSUM + PDIF*SDIF) C C------- dPsi/dm DZDM(JM) = DZDM(JM) + QOPI*(-PSUM*DSIM + PDIF*DSIM) DZDM(JO) = DZDM(JO) + QOPI*(-PSUM*DSIO - PDIF*DSIO) DZDM(JP) = DZDM(JP) + QOPI*( PSUM*(DSIO+DSIM) & + PDIF*(DSIO-DSIM)) C C------- dPsi/dni PSNI = PSX1*X1I + PSX0*(X1I+X2I)*0.5 + PSYY*YYI PDNI = PDX1*X1I + PDX0*(X1I+X2I)*0.5 + PDYY*YYI PSI_NI = PSI_NI + QOPI*(PSNI*SSUM + PDNI*SDIF) C QTANM = QTANM + QOPI*(PSNI*SSUM + PDNI*SDIF) C DQDM(JM) = DQDM(JM) + QOPI*(-PSNI*DSIM + PDNI*DSIM) DQDM(JO) = DQDM(JO) + QOPI*(-PSNI*DSIO - PDNI*DSIO) DQDM(JP) = DQDM(JP) + QOPI*( PSNI*(DSIO+DSIM) & + PDNI*(DSIO-DSIM)) C C C------- calculate source contribution to Psi for 0-2 half-panel DXINV = 1.0/(X0-X2) PSUM = X2*(T2-APAN) - X0*(T0-APAN) + 0.5*YY*(G0-G2) PDIF = ((X0+X2)*PSUM + RS0*(T0-APAN) - RS2*(T2-APAN) & + (X2-X0)*YY) * DXINV C PSX0 = -(T0-APAN) PSX2 = T2-APAN PSYY = 0.5*(G0-G2) C PDX0 = ((X0+X2)*PSX0 + PSUM + 2.0*X0*(T0-APAN) - PDIF) * DXINV PDX2 = ((X0+X2)*PSX2 + PSUM - 2.0*X2*(T2-APAN) + PDIF) * DXINV PDYY = ((X0+X2)*PSYY + 2.0*(X2-X0 + YY*(T0-T2)) ) * DXINV C DSP = SQRT((X(JQ)-X(JO))**2 + (Y(JQ)-Y(JO))**2) DSIP = 1.0/DSP C CCC SIG2 = (SIG(JQ) - SIG(JO))*DSIP CCC SIG0 = (SIG(JP) - SIG(JO))*DSIO CCC SSUM = SIG2 + SIG0 CCC SDIF = SIG2 - SIG0 C SSUM = (SIG(JQ) - SIG(JO))*DSIP + (SIG(JP) - SIG(JO))*DSIO SDIF = (SIG(JQ) - SIG(JO))*DSIP - (SIG(JP) - SIG(JO))*DSIO C PSI = PSI + QOPI*(PSUM*SSUM + PDIF*SDIF) C C------- dPsi/dm DZDM(JO) = DZDM(JO) + QOPI*(-PSUM*(DSIP+DSIO) & - PDIF*(DSIP-DSIO)) DZDM(JP) = DZDM(JP) + QOPI*( PSUM*DSIO - PDIF*DSIO) DZDM(JQ) = DZDM(JQ) + QOPI*( PSUM*DSIP + PDIF*DSIP) C C------- dPsi/dni PSNI = PSX0*(X1I+X2I)*0.5 + PSX2*X2I + PSYY*YYI PDNI = PDX0*(X1I+X2I)*0.5 + PDX2*X2I + PDYY*YYI PSI_NI = PSI_NI + QOPI*(PSNI*SSUM + PDNI*SDIF) C QTANM = QTANM + QOPI*(PSNI*SSUM + PDNI*SDIF) C DQDM(JO) = DQDM(JO) + QOPI*(-PSNI*(DSIP+DSIO) & - PDNI*(DSIP-DSIO)) DQDM(JP) = DQDM(JP) + QOPI*( PSNI*DSIO - PDNI*DSIO) DQDM(JQ) = DQDM(JQ) + QOPI*( PSNI*DSIP + PDNI*DSIP) C ENDIF C C------ calculate vortex panel contribution to Psi DXINV = 1.0/(X1-X2) PSIS = 0.5*X1*G1 - 0.5*X2*G2 + X2 - X1 + YY*(T1-T2) PSID = ((X1+X2)*PSIS + 0.5*(RS2*G2-RS1*G1 + X1*X1-X2*X2))*DXINV C PSX1 = 0.5*G1 PSX2 = -.5*G2 PSYY = T1-T2 C PDX1 = ((X1+X2)*PSX1 + PSIS - X1*G1 - PSID)*DXINV PDX2 = ((X1+X2)*PSX2 + PSIS + X2*G2 + PSID)*DXINV PDYY = ((X1+X2)*PSYY - YY*(G1-G2) )*DXINV C GSUM1 = GAMU(JP,1) + GAMU(JO,1) GSUM2 = GAMU(JP,2) + GAMU(JO,2) GDIF1 = GAMU(JP,1) - GAMU(JO,1) GDIF2 = GAMU(JP,2) - GAMU(JO,2) C GSUM = GAM(JP) + GAM(JO) GDIF = GAM(JP) - GAM(JO) C PSI = PSI + QOPI*(PSIS*GSUM + PSID*GDIF) C C------ dPsi/dGam DZDG(JO) = DZDG(JO) + QOPI*(PSIS-PSID) DZDG(JP) = DZDG(JP) + QOPI*(PSIS+PSID) C C------ dPsi/dni PSNI = PSX1*X1I + PSX2*X2I + PSYY*YYI PDNI = PDX1*X1I + PDX2*X2I + PDYY*YYI PSI_NI = PSI_NI + QOPI*(GSUM*PSNI + GDIF*PDNI) C QTAN1 = QTAN1 + QOPI*(GSUM1*PSNI + GDIF1*PDNI) QTAN2 = QTAN2 + QOPI*(GSUM2*PSNI + GDIF2*PDNI) C DQDG(JO) = DQDG(JO) + QOPI*(PSNI - PDNI) DQDG(JP) = DQDG(JP) + QOPI*(PSNI + PDNI) C IF(GEOLIN) THEN C C------- dPsi/dn DZDN(JO) = DZDN(JO)+ QOPI*GSUM*(PSX1*X1O + PSX2*X2O + PSYY*YYO) & + QOPI*GDIF*(PDX1*X1O + PDX2*X2O + PDYY*YYO) DZDN(JP) = DZDN(JP)+ QOPI*GSUM*(PSX1*X1P + PSX2*X2P + PSYY*YYP) & + QOPI*GDIF*(PDX1*X1P + PDX2*X2P + PDYY*YYP) C------- dPsi/dP Z_QDOF0 = Z_QDOF0 & + QOPI*((PSIS-PSID)*QF0(JO) + (PSIS+PSID)*QF0(JP)) Z_QDOF1 = Z_QDOF1 & + QOPI*((PSIS-PSID)*QF1(JO) + (PSIS+PSID)*QF1(JP)) Z_QDOF2 = Z_QDOF2 & + QOPI*((PSIS-PSID)*QF2(JO) + (PSIS+PSID)*QF2(JP)) Z_QDOF3 = Z_QDOF3 & + QOPI*((PSIS-PSID)*QF3(JO) + (PSIS+PSID)*QF3(JP)) ENDIF C C 10 CONTINUE C 11 CONTINUE PSIG = 0.5*YY*(G1-G2) + X2*(T2-APAN) - X1*(T1-APAN) PGAM = 0.5*X1*G1 - 0.5*X2*G2 + X2 - X1 + YY*(T1-T2) C PSIGX1 = -(T1-APAN) PSIGX2 = T2-APAN PSIGYY = 0.5*(G1-G2) PGAMX1 = 0.5*G1 PGAMX2 = -.5*G2 PGAMYY = T1-T2 C PSIGNI = PSIGX1*X1I + PSIGX2*X2I + PSIGYY*YYI PGAMNI = PGAMX1*X1I + PGAMX2*X2I + PGAMYY*YYI C C---- TE panel source and vortex strengths SIGTE1 = 0.5*SCS*(GAMU(JP,1) - GAMU(JO,1)) SIGTE2 = 0.5*SCS*(GAMU(JP,2) - GAMU(JO,2)) GAMTE1 = -.5*SDS*(GAMU(JP,1) - GAMU(JO,1)) GAMTE2 = -.5*SDS*(GAMU(JP,2) - GAMU(JO,2)) C SIGTE = 0.5*SCS*(GAM(JP) - GAM(JO)) GAMTE = -.5*SDS*(GAM(JP) - GAM(JO)) C C---- TE panel contribution to Psi PSI = PSI + HOPI*(PSIG*SIGTE + PGAM*GAMTE) C C---- dPsi/dGam DZDG(JO) = DZDG(JO) - HOPI*PSIG*SCS*0.5 DZDG(JP) = DZDG(JP) + HOPI*PSIG*SCS*0.5 C DZDG(JO) = DZDG(JO) + HOPI*PGAM*SDS*0.5 DZDG(JP) = DZDG(JP) - HOPI*PGAM*SDS*0.5 C C---- dPsi/dni PSI_NI = PSI_NI + HOPI*(PSIGNI*SIGTE + PGAMNI*GAMTE) C QTAN1 = QTAN1 + HOPI*(PSIGNI*SIGTE1 + PGAMNI*GAMTE1) QTAN2 = QTAN2 + HOPI*(PSIGNI*SIGTE2 + PGAMNI*GAMTE2) C DQDG(JO) = DQDG(JO) - HOPI*(PSIGNI*0.5*SCS - PGAMNI*0.5*SDS) DQDG(JP) = DQDG(JP) + HOPI*(PSIGNI*0.5*SCS - PGAMNI*0.5*SDS) C IF(GEOLIN) THEN C C----- dPsi/dn DZDN(JO) = DZDN(JO) & + HOPI*(PSIGX1*X1O + PSIGX2*X2O + PSIGYY*YYO)*SIGTE & + HOPI*(PGAMX1*X1O + PGAMX2*X2O + PGAMYY*YYO)*GAMTE DZDN(JP) = DZDN(JP) & + HOPI*(PSIGX1*X1P + PSIGX2*X2P + PSIGYY*YYP)*SIGTE & + HOPI*(PGAMX1*X1P + PGAMX2*X2P + PGAMYY*YYP)*GAMTE C C----- dPsi/dP Z_QDOF0 = Z_QDOF0 + HOPI*PSIG*0.5*(QF0(JP)-QF0(JO))*SCS & - HOPI*PGAM*0.5*(QF0(JP)-QF0(JO))*SDS Z_QDOF1 = Z_QDOF1 + HOPI*PSIG*0.5*(QF1(JP)-QF1(JO))*SCS & - HOPI*PGAM*0.5*(QF1(JP)-QF1(JO))*SDS Z_QDOF2 = Z_QDOF2 + HOPI*PSIG*0.5*(QF2(JP)-QF2(JO))*SCS & - HOPI*PGAM*0.5*(QF2(JP)-QF2(JO))*SDS Z_QDOF3 = Z_QDOF3 + HOPI*PSIG*0.5*(QF3(JP)-QF3(JO))*SCS & - HOPI*PGAM*0.5*(QF3(JP)-QF3(JO))*SDS C ENDIF C 12 CONTINUE C C**** Freestream terms PSI = PSI + QINF*(COSA*YI - SINA*XI) C C---- dPsi/dn PSI_NI = PSI_NI + QINF*(COSA*NYI - SINA*NXI) C QTAN1 = QTAN1 + QINF*NYI QTAN2 = QTAN2 - QINF*NXI C C---- dPsi/dQinf Z_QINF = Z_QINF + (COSA*YI - SINA*XI) C C---- dPsi/dalfa Z_ALFA = Z_ALFA - QINF*(SINA*YI + COSA*XI) C IF(.NOT.LIMAGE) RETURN C C C THETJ(0) = 0.0 DO 20 JO=1, N JP = JO+1 C JM = JO-1 JQ = JP+1 C IF(JO.EQ.1) THEN JM = JO ELSE IF(JO.EQ.N-1) THEN JQ = JP ELSE IF(JO.EQ.N) THEN JP = 1 IF((X(JO)-X(JP))**2 + (Y(JO)-Y(JP))**2 .LT. SEPS**2) GO TO 22 ENDIF C DSO = SQRT((X(JO)-X(JP))**2 + (Y(JO)-Y(JP))**2) C C------ skip null panel IF(DSO .EQ. 0.0) GO TO 20 C DSIO = 1.0 / DSO C ccc APAN = APANEL(JO) APAN = PI - APANEL(JO) + 2.0*ALFA C XJO = X(JO) + 2.0*(YIMAGE+Y(JO))*SINA YJO = Y(JO) - 2.0*(YIMAGE+Y(JO))*COSA XJP = X(JP) + 2.0*(YIMAGE+Y(JP))*SINA YJP = Y(JP) - 2.0*(YIMAGE+Y(JP))*COSA C RX1 = XI - XJO RY1 = YI - YJO RX2 = XI - XJP RY2 = YI - YJP C SX = (XJP - XJO) * DSIO SY = (YJP - YJO) * DSIO C X1 = SX*RX1 + SY*RY1 X2 = SX*RX2 + SY*RY2 YY = SX*RY1 - SY*RX1 C RS1 = RX1*RX1 + RY1*RY1 RS2 = RX2*RX2 + RY2*RY2 C C------ set log(r^2) and arctan(x/y), correcting for reflection if any G1 = LOG(RS1) ccc T1 = ATAN2(X1,YY) T1 = ATANC(X1,YY,THETJ(JO-1)) THETJ(JO) = T1 C G2 = LOG(RS2) ccc T2 = ATAN2(X2,YY) T2 = ATANC(X2,YY,THETJ(JO)) THETJ(JP) = T2 C X1I = SX*NXI + SY*NYI X2I = SX*NXI + SY*NYI YYI = SX*NYI - SY*NXI C IF(GEOLIN) THEN NXO = NX(JO) NYO = NY(JO) NXP = NX(JP) NYP = NY(JP) C X1O =-((RX1-X1*SX)*NXO + (RY1-X1*SY)*NYO)*DSIO-(SX*NXO+SY*NYO) X1P = ((RX1-X1*SX)*NXP + (RY1-X1*SY)*NYP)*DSIO X2O =-((RX2-X2*SX)*NXO + (RY2-X2*SY)*NYO)*DSIO X2P = ((RX2-X2*SX)*NXP + (RY2-X2*SY)*NYP)*DSIO-(SX*NXP+SY*NYP) YYO = ((RX1+X1*SY)*NYO - (RY1-X1*SX)*NXO)*DSIO-(SX*NYO-SY*NXO) YYP =-((RX1-X1*SY)*NYP - (RY1+X1*SX)*NXP)*DSIO ENDIF C IF(JO.EQ.N) GO TO 21 C IF(SIGLIN) THEN C C------- set up midpoint quantities X0 = 0.5*(X1+X2) RS0 = X0*X0 + YY*YY G0 = LOG(RS0) ccc T0 = ATAN2(X0,YY) T0 = ATANC(X0,YY,THETJ(JO)) C C------- calculate source contribution to Psi for 1-0 half-panel DXINV = 1.0/(X1-X0) PSUM = X0*(T0-APAN) - X1*(T1-APAN) + 0.5*YY*(G1-G0) PDIF = ((X1+X0)*PSUM + RS1*(T1-APAN) - RS0*(T0-APAN) & + (X0-X1)*YY) * DXINV C PSX1 = -(T1-APAN) PSX0 = T0-APAN PSYY = 0.5*(G1-G0) C PDX1 = ((X1+X0)*PSX1 + PSUM + 2.0*X1*(T1-APAN) - PDIF) * DXINV PDX0 = ((X1+X0)*PSX0 + PSUM - 2.0*X0*(T0-APAN) + PDIF) * DXINV PDYY = ((X1+X0)*PSYY + 2.0*(X0-X1 + YY*(T1-T0)) ) * DXINV C DSM = SQRT((X(JP)-X(JM))**2 + (Y(JP)-Y(JM))**2) DSIM = 1.0/DSM C CCC SIG0 = (SIG(JP) - SIG(JO))*DSIO CCC SIG1 = (SIG(JP) - SIG(JM))*DSIM CCC SSUM = SIG0 + SIG1 CCC SDIF = SIG0 - SIG1 C SSUM = (SIG(JP) - SIG(JO))*DSIO + (SIG(JP) - SIG(JM))*DSIM SDIF = (SIG(JP) - SIG(JO))*DSIO - (SIG(JP) - SIG(JM))*DSIM C PSI = PSI + QOPI*(PSUM*SSUM + PDIF*SDIF) C C------- dPsi/dm DZDM(JM) = DZDM(JM) + QOPI*(-PSUM*DSIM + PDIF*DSIM) DZDM(JO) = DZDM(JO) + QOPI*(-PSUM*DSIO - PDIF*DSIO) DZDM(JP) = DZDM(JP) + QOPI*( PSUM*(DSIO+DSIM) & + PDIF*(DSIO-DSIM)) C C------- dPsi/dni PSNI = PSX1*X1I + PSX0*(X1I+X2I)*0.5 + PSYY*YYI PDNI = PDX1*X1I + PDX0*(X1I+X2I)*0.5 + PDYY*YYI PSI_NI = PSI_NI + QOPI*(PSNI*SSUM + PDNI*SDIF) C QTANM = QTANM + QOPI*(PSNI*SSUM + PDNI*SDIF) C DQDM(JM) = DQDM(JM) + QOPI*(-PSNI*DSIM + PDNI*DSIM) DQDM(JO) = DQDM(JO) + QOPI*(-PSNI*DSIO - PDNI*DSIO) DQDM(JP) = DQDM(JP) + QOPI*( PSNI*(DSIO+DSIM) & + PDNI*(DSIO-DSIM)) C C C------- calculate source contribution to Psi for 0-2 half-panel DXINV = 1.0/(X0-X2) PSUM = X2*(T2-APAN) - X0*(T0-APAN) + 0.5*YY*(G0-G2) PDIF = ((X0+X2)*PSUM + RS0*(T0-APAN) - RS2*(T2-APAN) & + (X2-X0)*YY) * DXINV C PSX0 = -(T0-APAN) PSX2 = T2-APAN PSYY = 0.5*(G0-G2) C PDX0 = ((X0+X2)*PSX0 + PSUM + 2.0*X0*(T0-APAN) - PDIF) * DXINV PDX2 = ((X0+X2)*PSX2 + PSUM - 2.0*X2*(T2-APAN) + PDIF) * DXINV PDYY = ((X0+X2)*PSYY + 2.0*(X2-X0 + YY*(T0-T2)) ) * DXINV C DSP = SQRT((X(JQ)-X(JO))**2 + (Y(JQ)-Y(JO))**2) DSIP = 1.0/DSP C CCC SIG2 = (SIG(JQ) - SIG(JO))*DSIP CCC SIG0 = (SIG(JP) - SIG(JO))*DSIO CCC SSUM = SIG2 + SIG0 CCC SDIF = SIG2 - SIG0 C SSUM = (SIG(JQ) - SIG(JO))*DSIP + (SIG(JP) - SIG(JO))*DSIO SDIF = (SIG(JQ) - SIG(JO))*DSIP - (SIG(JP) - SIG(JO))*DSIO C PSI = PSI + QOPI*(PSUM*SSUM + PDIF*SDIF) C C------- dPsi/dm DZDM(JO) = DZDM(JO) + QOPI*(-PSUM*(DSIP+DSIO) & - PDIF*(DSIP-DSIO)) DZDM(JP) = DZDM(JP) + QOPI*( PSUM*DSIO - PDIF*DSIO) DZDM(JQ) = DZDM(JQ) + QOPI*( PSUM*DSIP + PDIF*DSIP) C C------- dPsi/dni PSNI = PSX0*(X1I+X2I)*0.5 + PSX2*X2I + PSYY*YYI PDNI = PDX0*(X1I+X2I)*0.5 + PDX2*X2I + PDYY*YYI PSI_NI = PSI_NI + QOPI*(PSNI*SSUM + PDNI*SDIF) C QTANM = QTANM + QOPI*(PSNI*SSUM + PDNI*SDIF) C DQDM(JO) = DQDM(JO) + QOPI*(-PSNI*(DSIP+DSIO) & - PDNI*(DSIP-DSIO)) DQDM(JP) = DQDM(JP) + QOPI*( PSNI*DSIO - PDNI*DSIO) DQDM(JQ) = DQDM(JQ) + QOPI*( PSNI*DSIP + PDNI*DSIP) C ENDIF C C------ calculate vortex panel contribution to Psi DXINV = 1.0/(X1-X2) PSIS = 0.5*X1*G1 - 0.5*X2*G2 + X2 - X1 + YY*(T1-T2) PSID = ((X1+X2)*PSIS + 0.5*(RS2*G2-RS1*G1 + X1*X1-X2*X2))*DXINV C PSX1 = 0.5*G1 PSX2 = -.5*G2 PSYY = T1-T2 C PDX1 = ((X1+X2)*PSX1 + PSIS - X1*G1 - PSID)*DXINV PDX2 = ((X1+X2)*PSX2 + PSIS + X2*G2 + PSID)*DXINV PDYY = ((X1+X2)*PSYY - YY*(G1-G2) )*DXINV C GSUM1 = GAMU(JP,1) + GAMU(JO,1) GSUM2 = GAMU(JP,2) + GAMU(JO,2) GDIF1 = GAMU(JP,1) - GAMU(JO,1) GDIF2 = GAMU(JP,2) - GAMU(JO,2) C GSUM = GAM(JP) + GAM(JO) GDIF = GAM(JP) - GAM(JO) C PSI = PSI - QOPI*(PSIS*GSUM + PSID*GDIF) C C------ dPsi/dGam DZDG(JO) = DZDG(JO) - QOPI*(PSIS-PSID) DZDG(JP) = DZDG(JP) - QOPI*(PSIS+PSID) C C------ dPsi/dni PSNI = PSX1*X1I + PSX2*X2I + PSYY*YYI PDNI = PDX1*X1I + PDX2*X2I + PDYY*YYI PSI_NI = PSI_NI - QOPI*(GSUM*PSNI + GDIF*PDNI) C QTAN1 = QTAN1 - QOPI*(GSUM1*PSNI + GDIF1*PDNI) QTAN2 = QTAN2 - QOPI*(GSUM2*PSNI + GDIF2*PDNI) C DQDG(JO) = DQDG(JO) - QOPI*(PSNI - PDNI) DQDG(JP) = DQDG(JP) - QOPI*(PSNI + PDNI) C IF(GEOLIN) THEN C C------- dPsi/dn DZDN(JO) = DZDN(JO)- QOPI*GSUM*(PSX1*X1O + PSX2*X2O + PSYY*YYO) & - QOPI*GDIF*(PDX1*X1O + PDX2*X2O + PDYY*YYO) DZDN(JP) = DZDN(JP)- QOPI*GSUM*(PSX1*X1P + PSX2*X2P + PSYY*YYP) & - QOPI*GDIF*(PDX1*X1P + PDX2*X2P + PDYY*YYP) C------- dPsi/dP Z_QDOF0 = Z_QDOF0 & - QOPI*((PSIS-PSID)*QF0(JO) + (PSIS+PSID)*QF0(JP)) Z_QDOF1 = Z_QDOF1 & - QOPI*((PSIS-PSID)*QF1(JO) + (PSIS+PSID)*QF1(JP)) Z_QDOF2 = Z_QDOF2 & - QOPI*((PSIS-PSID)*QF2(JO) + (PSIS+PSID)*QF2(JP)) Z_QDOF3 = Z_QDOF3 & - QOPI*((PSIS-PSID)*QF3(JO) + (PSIS+PSID)*QF3(JP)) ENDIF C C 20 CONTINUE C 21 CONTINUE PSIG = 0.5*YY*(G1-G2) + X2*(T2-APAN) - X1*(T1-APAN) PGAM = 0.5*X1*G1 - 0.5*X2*G2 + X2 - X1 + YY*(T1-T2) C PSIGX1 = -(T1-APAN) PSIGX2 = T2-APAN PSIGYY = 0.5*(G1-G2) PGAMX1 = 0.5*G1 PGAMX2 = -.5*G2 PGAMYY = T1-T2 C PSIGNI = PSIGX1*X1I + PSIGX2*X2I + PSIGYY*YYI PGAMNI = PGAMX1*X1I + PGAMX2*X2I + PGAMYY*YYI C C---- TE panel source and vortex strengths SIGTE1 = 0.5*SCS*(GAMU(JP,1) - GAMU(JO,1)) SIGTE2 = 0.5*SCS*(GAMU(JP,2) - GAMU(JO,2)) GAMTE1 = -.5*SDS*(GAMU(JP,1) - GAMU(JO,1)) GAMTE2 = -.5*SDS*(GAMU(JP,2) - GAMU(JO,2)) C SIGTE = 0.5*SCS*(GAM(JP) - GAM(JO)) GAMTE = -.5*SDS*(GAM(JP) - GAM(JO)) C C---- TE panel contribution to Psi PSI = PSI + HOPI*(PSIG*SIGTE - PGAM*GAMTE) C C---- dPsi/dGam DZDG(JO) = DZDG(JO) - HOPI*PSIG*SCS*0.5 DZDG(JP) = DZDG(JP) + HOPI*PSIG*SCS*0.5 C DZDG(JO) = DZDG(JO) - HOPI*PGAM*SDS*0.5 DZDG(JP) = DZDG(JP) + HOPI*PGAM*SDS*0.5 C C---- dPsi/dni PSI_NI = PSI_NI + HOPI*(PSIGNI*SIGTE - PGAMNI*GAMTE) C QTAN1 = QTAN1 + HOPI*(PSIGNI*SIGTE1 - PGAMNI*GAMTE1) QTAN2 = QTAN2 + HOPI*(PSIGNI*SIGTE2 - PGAMNI*GAMTE2) C DQDG(JO) = DQDG(JO) - HOPI*(PSIGNI*0.5*SCS + PGAMNI*0.5*SDS) DQDG(JP) = DQDG(JP) + HOPI*(PSIGNI*0.5*SCS + PGAMNI*0.5*SDS) C IF(GEOLIN) THEN C C----- dPsi/dn DZDN(JO) = DZDN(JO) & + HOPI*(PSIGX1*X1O + PSIGX2*X2O + PSIGYY*YYO)*SIGTE & - HOPI*(PGAMX1*X1O + PGAMX2*X2O + PGAMYY*YYO)*GAMTE DZDN(JP) = DZDN(JP) & + HOPI*(PSIGX1*X1P + PSIGX2*X2P + PSIGYY*YYP)*SIGTE & - HOPI*(PGAMX1*X1P + PGAMX2*X2P + PGAMYY*YYP)*GAMTE C C----- dPsi/dP Z_QDOF0 = Z_QDOF0 + HOPI*PSIG*0.5*(QF0(JP)-QF0(JO))*SCS & + HOPI*PGAM*0.5*(QF0(JP)-QF0(JO))*SDS Z_QDOF1 = Z_QDOF1 + HOPI*PSIG*0.5*(QF1(JP)-QF1(JO))*SCS & + HOPI*PGAM*0.5*(QF1(JP)-QF1(JO))*SDS Z_QDOF2 = Z_QDOF2 + HOPI*PSIG*0.5*(QF2(JP)-QF2(JO))*SCS & + HOPI*PGAM*0.5*(QF2(JP)-QF2(JO))*SDS Z_QDOF3 = Z_QDOF3 + HOPI*PSIG*0.5*(QF3(JP)-QF3(JO))*SCS & + HOPI*PGAM*0.5*(QF3(JP)-QF3(JO))*SDS C ENDIF C 22 CONTINUE C RETURN END SUBROUTINE PSWLIN(I,XI,YI,NXI,NYI,PSI,PSI_NI) C-------------------------------------------------------------------- C Calculates current streamfunction Psi and tangential velocity C Qtan at panel node or wake node I due to freestream and wake C sources Sig. Also calculates sensitivity vectors dPsi/dSig C (DZDM) and dQtan/dSig (DQDM). C C Airfoil: 1 < I < N C Wake: N+1 < I < N+NW C-------------------------------------------------------------------- INCLUDE 'XFOIL.INC' REAL NXI, NYI REAL THETJ(0:IZX) C IO = I C COSA = COS(ALFA) SINA = SIN(ALFA) C DO 4 JO=N+1, N+NW DZDM(JO) = 0.0 DQDM(JO) = 0.0 4 CONTINUE C PSI = 0. PSI_NI = 0. C THETJ(0) = 0.0 DO 20 JO=N+1, N+NW-1 C JP = JO+1 C JM = JO-1 JQ = JP+1 IF(JO.EQ.N+1) THEN JM = JO ELSE IF(JO.EQ.N+NW-1) THEN JQ = JP ENDIF C DSO = SQRT((X(JO)-X(JP))**2 + (Y(JO)-Y(JP))**2) DSIO = 1.0 / DSO C APAN = APANEL(JO) C RX1 = XI - X(JO) RY1 = YI - Y(JO) RX2 = XI - X(JP) RY2 = YI - Y(JP) C SX = (X(JP) - X(JO)) * DSIO SY = (Y(JP) - Y(JO)) * DSIO C X1 = SX*RX1 + SY*RY1 X2 = SX*RX2 + SY*RY2 YY = SX*RY1 - SY*RX1 C RS1 = RX1*RX1 + RY1*RY1 RS2 = RX2*RX2 + RY2*RY2 C IF(IO.NE.JO .AND. RS1.GT.0.0) THEN G1 = LOG(RS1) ccc T1 = ATAN2(X1,YY) T1 = ATANC(X1,YY,THETJ(JO-1)) ELSE G1 = 0.0 T1 = 0.0 ENDIF THETJ(JO) = T1 C IF(IO.NE.JP .AND. RS2.GT.0.0) THEN G2 = LOG(RS2) ccc T2 = ATAN2(X2,YY) T2 = ATANC(X2,YY,THETJ(JO)) ELSE G2 = 0.0 T2 = 0.0 ENDIF THETJ(JP) = T2 C X1I = SX*NXI + SY*NYI X2I = SX*NXI + SY*NYI YYI = SX*NYI - SY*NXI C C------- set up midpoint quantities X0 = 0.5*(X1+X2) RS0 = X0*X0 + YY*YY G0 = LOG(RS0) ccc T0 = ATAN2(X0,YY) T0 = ATANC(X0,YY,THETJ(JO-1)) C C------- calculate source contribution to Psi for 1-0 half-panel DXINV = 1.0/(X1-X0) PSUM = X0*(T0-APAN) - X1*(T1-APAN) + 0.5*YY*(G1-G0) PDIF = ((X1+X0)*PSUM + RS1*(T1-APAN) - RS0*(T0-APAN) & + (X0-X1)*YY) * DXINV C PSX1 = -(T1-APAN) PSX0 = T0-APAN PSYY = 0.5*(G1-G0) C PDX1 = ((X1+X0)*PSX1 + PSUM + 2.0*X1*(T1-APAN) - PDIF) * DXINV PDX0 = ((X1+X0)*PSX0 + PSUM - 2.0*X0*(T0-APAN) + PDIF) * DXINV PDYY = ((X1+X0)*PSYY + 2.0*(X0-X1 + YY*(T1-T0)) ) * DXINV C DSM = SQRT((X(JP)-X(JM))**2 + (Y(JP)-Y(JM))**2) DSIM = 1.0/DSM C CCC SIG0 = (SIG(JP) - SIG(JO))*DSIO CCC SIG1 = (SIG(JP) - SIG(JM))*DSIM CCC SSUM = SIG0 + SIG1 CCC SDIF = SIG0 - SIG1 C SSUM = (SIG(JP) - SIG(JO))*DSIO + (SIG(JP) - SIG(JM))*DSIM SDIF = (SIG(JP) - SIG(JO))*DSIO - (SIG(JP) - SIG(JM))*DSIM C PSI = PSI + QOPI*(PSUM*SSUM + PDIF*SDIF) C C------- dPsi/dm DZDM(JM) = DZDM(JM) + QOPI*(-PSUM*DSIM + PDIF*DSIM) DZDM(JO) = DZDM(JO) + QOPI*(-PSUM*DSIO - PDIF*DSIO) DZDM(JP) = DZDM(JP) + QOPI*( PSUM*(DSIO+DSIM) & + PDIF*(DSIO-DSIM)) C C------- dPsi/dni PSNI = PSX1*X1I + PSX0*(X1I+X2I)*0.5 + PSYY*YYI PDNI = PDX1*X1I + PDX0*(X1I+X2I)*0.5 + PDYY*YYI PSI_NI = PSI_NI + QOPI*(PSNI*SSUM + PDNI*SDIF) C DQDM(JM) = DQDM(JM) + QOPI*(-PSNI*DSIM + PDNI*DSIM) DQDM(JO) = DQDM(JO) + QOPI*(-PSNI*DSIO - PDNI*DSIO) DQDM(JP) = DQDM(JP) + QOPI*( PSNI*(DSIO+DSIM) & + PDNI*(DSIO-DSIM)) C C C------- calculate source contribution to Psi for 0-2 half-panel DXINV = 1.0/(X0-X2) PSUM = X2*(T2-APAN) - X0*(T0-APAN) + 0.5*YY*(G0-G2) PDIF = ((X0+X2)*PSUM + RS0*(T0-APAN) - RS2*(T2-APAN) & + (X2-X0)*YY) * DXINV C PSX0 = -(T0-APAN) PSX2 = T2-APAN PSYY = 0.5*(G0-G2) C PDX0 = ((X0+X2)*PSX0 + PSUM + 2.0*X0*(T0-APAN) - PDIF) * DXINV PDX2 = ((X0+X2)*PSX2 + PSUM - 2.0*X2*(T2-APAN) + PDIF) * DXINV PDYY = ((X0+X2)*PSYY + 2.0*(X2-X0 + YY*(T0-T2)) ) * DXINV C DSP = SQRT((X(JQ)-X(JO))**2 + (Y(JQ)-Y(JO))**2) DSIP = 1.0/DSP C CCC SIG2 = (SIG(JQ) - SIG(JO))*DSIP CCC SIG0 = (SIG(JP) - SIG(JO))*DSIO CCC SSUM = SIG2 + SIG0 CCC SDIF = SIG2 - SIG0 C SSUM = (SIG(JQ) - SIG(JO))*DSIP + (SIG(JP) - SIG(JO))*DSIO SDIF = (SIG(JQ) - SIG(JO))*DSIP - (SIG(JP) - SIG(JO))*DSIO C PSI = PSI + QOPI*(PSUM*SSUM + PDIF*SDIF) C C------- dPsi/dm DZDM(JO) = DZDM(JO) + QOPI*(-PSUM*(DSIP+DSIO) & - PDIF*(DSIP-DSIO)) DZDM(JP) = DZDM(JP) + QOPI*( PSUM*DSIO - PDIF*DSIO) DZDM(JQ) = DZDM(JQ) + QOPI*( PSUM*DSIP + PDIF*DSIP) C C------- dPsi/dni PSNI = PSX0*(X1I+X2I)*0.5 + PSX2*X2I + PSYY*YYI PDNI = PDX0*(X1I+X2I)*0.5 + PDX2*X2I + PDYY*YYI PSI_NI = PSI_NI + QOPI*(PSNI*SSUM + PDNI*SDIF) C DQDM(JO) = DQDM(JO) + QOPI*(-PSNI*(DSIP+DSIO) & - PDNI*(DSIP-DSIO)) DQDM(JP) = DQDM(JP) + QOPI*( PSNI*DSIO - PDNI*DSIO) DQDM(JQ) = DQDM(JQ) + QOPI*( PSNI*DSIP + PDNI*DSIP) C 20 CONTINUE C RETURN END SUBROUTINE GGCALC C-------------------------------------------------------------- C Calculates two surface vorticity (gamma) distributions C for alpha = 0, 90 degrees. These are superimposed C in SPECAL or SPECCL for specified alpha or CL. C-------------------------------------------------------------- INCLUDE 'XFOIL.INC' C C---- distance of internal control point ahead of sharp TE C- (fraction of smaller panel length adjacent to TE) BWT = 0.1 C WRITE(*,*) 'Calculating unit vorticity distributions ...' C DO 10 I=1, N GAM(I) = 0. GAMU(I,1) = 0. GAMU(I,2) = 0. 10 CONTINUE PSIO = 0. C C---- Set up matrix system for Psi = Psio on airfoil surface. C- The unknowns are (dGamma)i and dPsio. DO 20 I=1, N C C------ calculate Psi and dPsi/dGamma array for current node CALL PSILIN(I,X(I),Y(I),NX(I),NY(I),PSI,PSI_N,.FALSE.,.TRUE.) C PSIINF = QINF*(COS(ALFA)*Y(I) - SIN(ALFA)*X(I)) C C------ RES1 = PSI( 0) - PSIO C------ RES2 = PSI(90) - PSIO RES1 = QINF*Y(I) RES2 = -QINF*X(I) C C------ dRes/dGamma DO 201 J=1, N AIJ(I,J) = DZDG(J) 201 CONTINUE C DO 202 J=1, N BIJ(I,J) = -DZDM(J) 202 CONTINUE C C------ dRes/dPsio AIJ(I,N+1) = -1.0 C GAMU(I,1) = -RES1 GAMU(I,2) = -RES2 C 20 CONTINUE C C---- set Kutta condition C- RES = GAM(1) + GAM(N) RES = 0. C DO 30 J=1, N+1 AIJ(N+1,J) = 0.0 30 CONTINUE C AIJ(N+1,1) = 1.0 AIJ(N+1,N) = 1.0 C GAMU(N+1,1) = -RES GAMU(N+1,2) = -RES C C---- set up Kutta condition (no direct source influence) DO 32 J=1, N BIJ(N+1,J) = 0. 32 CONTINUE C IF(SHARP) THEN C----- set zero internal velocity in TE corner C C----- set TE bisector angle AG1 = ATAN2(-YP(1),-XP(1) ) AG2 = ATANC( YP(N), XP(N),AG1) ABIS = 0.5*(AG1+AG2) CBIS = COS(ABIS) SBIS = SIN(ABIS) C C----- minimum panel length adjacent to TE DS1 = SQRT( (X(1)-X(2) )**2 + (Y(1)-Y(2) )**2 ) DS2 = SQRT( (X(N)-X(N-1))**2 + (Y(N)-Y(N-1))**2 ) DSMIN = MIN( DS1 , DS2 ) C C----- control point on bisector just ahead of TE point XBIS = XTE - BWT*DSMIN*CBIS YBIS = YTE - BWT*DSMIN*SBIS ccc write(*,*) xbis, ybis C C----- set velocity component along bisector line CALL PSILIN(0,XBIS,YBIS,-SBIS,CBIS,PSI,QBIS,.FALSE.,.TRUE.) C CCC--- RES = DQDGj*Gammaj + DQDMj*Massj + QINF*(COSA*CBIS + SINA*SBIS) RES = QBIS C C----- dRes/dGamma DO J=1, N AIJ(N,J) = DQDG(J) ENDDO C C----- -dRes/dMass DO J=1, N BIJ(N,J) = -DQDM(J) ENDDO C C----- dRes/dPsio AIJ(N,N+1) = 0. C C----- -dRes/dUinf GAMU(N,1) = -CBIS C C----- -dRes/dVinf GAMU(N,2) = -SBIS C ENDIF C C---- LU-factor coefficient matrix AIJ CALL LUDCMP(IQX,N+1,AIJ,AIJPIV) LQAIJ = .TRUE. C C---- solve system for the two vorticity distributions CALL BAKSUB(IQX,N+1,AIJ,AIJPIV,GAMU(1,1)) CALL BAKSUB(IQX,N+1,AIJ,AIJPIV,GAMU(1,2)) C C---- set inviscid alpha=0,90 surface speeds for this geometry DO 50 I=1, N QINVU(I,1) = GAMU(I,1) QINVU(I,2) = GAMU(I,2) 50 CONTINUE C LGAMU = .TRUE. C RETURN END SUBROUTINE QWCALC C--------------------------------------------------------------- C Sets inviscid tangential velocity for alpha = 0, 90 C on wake due to freestream and airfoil surface vorticity. C--------------------------------------------------------------- INCLUDE 'XFOIL.INC' C C---- first wake point (same as TE) QINVU(N+1,1) = QINVU(N,1) QINVU(N+1,2) = QINVU(N,2) C C---- rest of wake DO 10 I=N+2, N+NW CALL PSILIN(I,X(I),Y(I),NX(I),NY(I),PSI,PSI_NI,.FALSE.,.FALSE.) QINVU(I,1) = QTAN1 QINVU(I,2) = QTAN2 10 CONTINUE C RETURN END SUBROUTINE QDCALC C----------------------------------------------------- C Calculates source panel influence coefficient C matrix for current airfoil and wake geometry. C----------------------------------------------------- INCLUDE 'XFOIL.INC' C WRITE(*,*) 'Calculating source influence matrix ...' C IF(.NOT.LADIJ) THEN C C----- calculate source influence matrix for airfoil surface if it doesn't exist DO 10 J=1, N C C------- multiply each dPsi/Sig vector by inverse of factored dPsi/dGam matrix CALL BAKSUB(IQX,N+1,AIJ,AIJPIV,BIJ(1,J)) C C------- store resulting dGam/dSig = dQtan/dSig vector DO 105 I=1, N DIJ(I,J) = BIJ(I,J) 105 CONTINUE C 10 CONTINUE LADIJ = .TRUE. C ENDIF C C---- set up coefficient matrix of dPsi/dm on airfoil surface DO 20 I=1, N CALL PSWLIN(I,X(I),Y(I),NX(I),NY(I),PSI,PSI_N) DO 202 J=N+1, N+NW BIJ(I,J) = -DZDM(J) 202 CONTINUE 20 CONTINUE C C---- set up Kutta condition (no direct source influence) DO 32 J=N+1, N+NW BIJ(N+1,J) = 0. 32 CONTINUE C C---- sharp TE gamma extrapolation also has no source influence IF(SHARP) THEN DO 34 J=N+1, N+NW BIJ(N,J) = 0. 34 CONTINUE ENDIF C C---- multiply by inverse of factored dPsi/dGam matrix DO 40 J=N+1, N+NW CALL BAKSUB(IQX,N+1,AIJ,AIJPIV,BIJ(1,J)) 40 CONTINUE C C---- set the source influence matrix for the wake sources DO 50 I=1, N DO 510 J=N+1, N+NW DIJ(I,J) = BIJ(I,J) 510 CONTINUE 50 CONTINUE C C**** Now we need to calculate the influence of sources on the wake velocities C C---- calculcate dQtan/dGam and dQtan/dSig at the wake points DO 70 I=N+1, N+NW C IW = I-N C C------ airfoil contribution at wake panel node CALL PSILIN(I,X(I),Y(I),NX(I),NY(I),PSI,PSI_N,.FALSE.,.TRUE.) C DO 710 J=1, N CIJ(IW,J) = DQDG(J) 710 CONTINUE C DO 720 J=1, N DIJ(I,J) = DQDM(J) 720 CONTINUE C C------ wake contribution CALL PSWLIN(I,X(I),Y(I),NX(I),NY(I),PSI,PSI_N) C DO 730 J=N+1, N+NW DIJ(I,J) = DQDM(J) 730 CONTINUE C 70 CONTINUE C C---- add on effect of all sources on airfoil vorticity which effects wake Qtan DO 80 I=N+1, N+NW IW = I-N C C------ airfoil surface source contribution first DO 810 J=1, N SUM = 0. DO 8100 K=1, N SUM = SUM + CIJ(IW,K)*DIJ(K,J) 8100 CONTINUE DIJ(I,J) = DIJ(I,J) + SUM 810 CONTINUE C C------ wake source contribution next DO 820 J=N+1, N+NW SUM = 0. DO 8200 K=1, N SUM = SUM + CIJ(IW,K)*BIJ(K,J) 8200 CONTINUE DIJ(I,J) = DIJ(I,J) + SUM 820 CONTINUE C 80 CONTINUE C C---- make sure first wake point has same velocity as trailing edge DO 90 J=1, N+NW DIJ(N+1,J) = DIJ(N,J) 90 CONTINUE C LWDIJ = .TRUE. C RETURN END SUBROUTINE XYWAKE C----------------------------------------------------- C Sets wake coordinate array for current surface C vorticity and/or mass source distributions. C----------------------------------------------------- INCLUDE 'XFOIL.INC' C WRITE(*,*) 'Calculating wake trajectory ...' C C---- number of wake points NW = N/8 + 2 IF(NW.GT.IWX) THEN WRITE(*,*) & 'Array size (IWX) too small. Last wake point index reduced.' NW = IWX ENDIF C DS1 = 0.5*(S(2) - S(1) + S(N) - S(N-1)) CALL SETEXP(SNEW(N+1),DS1,WAKLEN*CHORD,NW) C XTE = 0.5*(X(1)+X(N)) YTE = 0.5*(Y(1)+Y(N)) C C---- set first wake point a tiny distance behind TE I = N+1 SX = 0.5*(YP(N) - YP(1)) SY = 0.5*(XP(1) - XP(N)) SMOD = SQRT(SX**2 + SY**2) NX(I) = SX / SMOD NY(I) = SY / SMOD X(I) = XTE - 0.0001*NY(I) Y(I) = YTE + 0.0001*NX(I) S(I) = S(N) C C---- calculate streamfunction gradient components at first point CALL PSILIN(I,X(I),Y(I),1.0,0.0,PSI,PSI_X,.FALSE.,.FALSE.) CALL PSILIN(I,X(I),Y(I),0.0,1.0,PSI,PSI_Y,.FALSE.,.FALSE.) C C---- set unit vector normal to wake at first point NX(I+1) = -PSI_X / SQRT(PSI_X**2 + PSI_Y**2) NY(I+1) = -PSI_Y / SQRT(PSI_X**2 + PSI_Y**2) C C---- set angle of wake panel normal APANEL(I) = ATAN2( PSI_Y , PSI_X ) C C---- set rest of wake points DO 10 I=N+2, N+NW DS = SNEW(I) - SNEW(I-1) C C------ set new point DS downstream of last point X(I) = X(I-1) - DS*NY(I) Y(I) = Y(I-1) + DS*NX(I) S(I) = S(I-1) + DS C IF(I.EQ.N+NW) GO TO 10 C C------- calculate normal vector for next point CALL PSILIN(I,X(I),Y(I),1.0,0.0,PSI,PSI_X,.FALSE.,.FALSE.) CALL PSILIN(I,X(I),Y(I),0.0,1.0,PSI,PSI_Y,.FALSE.,.FALSE.) C NX(I+1) = -PSI_X / SQRT(PSI_X**2 + PSI_Y**2) NY(I+1) = -PSI_Y / SQRT(PSI_X**2 + PSI_Y**2) C C------- set angle of wake panel normal ccc APANEL(I) = ATAN2( PSI_Y , PSI_X ) APANEL(I) = ATANC( PSI_Y , PSI_X , APANEL(I-1) ) C 10 CONTINUE C C---- set wake presence flag and corresponding alpha LWAKE = .TRUE. AWAKE = ALFA C C---- old source influence matrix is invalid for the new wake geometry LWDIJ = .FALSE. C RETURN END SUBROUTINE STFIND C----------------------------------------- C Locates stagnation point arc length C location SST and panel index IST. C----------------------------------------- INCLUDE 'XFOIL.INC' C DO 10 I=1, N-1 IF(GAM(I).GE.0.0 .AND. GAM(I+1).LT.0.0) GO TO 11 10 CONTINUE C WRITE(*,*) 'STFIND: Stagnation point not found. Continuing ...' I = N/2 C 11 CONTINUE C IST = I DGAM = GAM(I+1) - GAM(I) DS = S(I+1) - S(I) C C---- evaluate so as to minimize roundoff for very small GAM(I) or GAM(I+1) IF(GAM(I) .LT. -GAM(I+1)) THEN SST = S(I) - DS*(GAM(I) /DGAM) ELSE SST = S(I+1) - DS*(GAM(I+1)/DGAM) ENDIF C C---- tweak stagnation point if it falls right on a node (very unlikely) IF(SST .LE. S(I) ) SST = S(I) + 1.0E-7 IF(SST .GE. S(I+1)) SST = S(I+1) - 1.0E-7 C SST_GO = (SST - S(I+1))/DGAM SST_GP = (S(I) - SST )/DGAM C RETURN END SUBROUTINE IBLPAN C------------------------------------------------------------- C Sets BL location -> panel location pointer array IPAN C------------------------------------------------------------- INCLUDE 'XFOIL.INC' C C---- top surface first IS = 1 C IBL = 1 DO 10 I=IST, 1, -1 IBL = IBL+1 IPAN(IBL,IS) = I VTI(IBL,IS) = 1.0 10 CONTINUE C IBLTE(IS) = IBL NBL(IS) = IBL C C---- bottom surface next IS = 2 C IBL = 1 DO 20 I=IST+1, N IBL = IBL+1 IPAN(IBL,IS) = I VTI(IBL,IS) = -1.0 20 CONTINUE C C---- wake IBLTE(IS) = IBL C DO 25 IW=1, NW I = N+IW IBL = IBLTE(IS)+IW IPAN(IBL,IS) = I VTI(IBL,IS) = -1.0 25 CONTINUE C NBL(IS) = IBLTE(IS) + NW C C---- upper wake pointers (for plotting only) DO 35 IW=1, NW IPAN(IBLTE(1)+IW,1) = IPAN(IBLTE(2)+IW,2) VTI(IBLTE(1)+IW,1) = 1.0 35 CONTINUE C C IBLMAX = MAX(IBLTE(1),IBLTE(2)) + NW IF(IBLMAX.GT.IVX) THEN WRITE(*,*) ' *** BL array overflow.' WRITE(*,*) ' *** Increase IVX to at least', IBLMAX STOP ENDIF C LIPAN = .TRUE. RETURN END SUBROUTINE XICALC C------------------------------------------------------------- C Sets BL arc length array on each airfoil side and wake C------------------------------------------------------------- INCLUDE 'XFOIL.INC' DATA XFEPS / 1.0E-7 / C C---- minimum xi node arc length near stagnation point XEPS = XFEPS*(S(N)-S(1)) C IS = 1 C XSSI(1,IS) = 0. C DO 10 IBL=2, IBLTE(IS) I = IPAN(IBL,IS) XSSI(IBL,IS) = MAX( SST - S(I) , XEPS ) 10 CONTINUE C C IS = 2 C XSSI(1,IS) = 0. C DO 20 IBL=2, IBLTE(IS) I = IPAN(IBL,IS) XSSI(IBL,IS) = MAX( S(I) - SST , XEPS ) 20 CONTINUE C C IS1 = 1 IS2 = 2 C IBL1 = IBLTE(IS1) + 1 XSSI(IBL1,IS1) = XSSI(IBL1-1,IS1) C IBL2 = IBLTE(IS2) + 1 XSSI(IBL2,IS2) = XSSI(IBL2-1,IS2) C DO 25 IBL=IBLTE(IS)+2, NBL(IS) I = IPAN(IBL,IS) DXSSI = SQRT((X(I)-X(I-1))**2 + (Y(I)-Y(I-1))**2) C IBL1 = IBLTE(IS1) + IBL - IBLTE(IS) IBL2 = IBLTE(IS2) + IBL - IBLTE(IS) XSSI(IBL1,IS1) = XSSI(IBL1-1,IS1) + DXSSI XSSI(IBL2,IS2) = XSSI(IBL2-1,IS2) + DXSSI 25 CONTINUE C C---- trailing edge flap length to TE gap ratio TELRAT = 2.50 C C---- set up parameters for TE flap cubics C ccc DWDXTE = YP(1)/XP(1) + YP(N)/XP(N) !!! BUG 2/2/95 C CROSP = (XP(1)*YP(N) - YP(1)*XP(N)) & / SQRT( (XP(1)**2 + YP(1)**2) & *(XP(N)**2 + YP(N)**2) ) DWDXTE = CROSP / SQRT(1.0 - CROSP**2) C C---- limit cubic to avoid absurd TE gap widths DWDXTE = MAX(DWDXTE,-3.0/TELRAT) DWDXTE = MIN(DWDXTE, 3.0/TELRAT) C AA = 3.0 + TELRAT*DWDXTE BB = -2.0 - TELRAT*DWDXTE C IF(SHARP) THEN DO 30 IW=1, NW WGAP(IW) = 0. 30 CONTINUE ELSE C----- set TE flap (wake gap) array IS = 2 DO 35 IW=1, NW IBL = IBLTE(IS) + IW ZN = 1.0 - (XSSI(IBL,IS)-XSSI(IBLTE(IS),IS)) / (TELRAT*ANTE) WGAP(IW) = 0. IF(ZN.GE.0.0) WGAP(IW) = ANTE * (AA + BB*ZN)*ZN**2 35 CONTINUE ENDIF C RETURN END SUBROUTINE UICALC C-------------------------------------------------------------- C Sets inviscid Ue from panel inviscid tangential velocity C-------------------------------------------------------------- INCLUDE 'XFOIL.INC' C DO 10 IS=1, 2 UINV (1,IS) = 0. UINV_A(1,IS) = 0. DO 110 IBL=2, NBL(IS) I = IPAN(IBL,IS) UINV (IBL,IS) = VTI(IBL,IS)*QINV (I) UINV_A(IBL,IS) = VTI(IBL,IS)*QINV_A(I) 110 CONTINUE 10 CONTINUE C RETURN END SUBROUTINE UECALC C-------------------------------------------------------------- C Sets viscous Ue from panel viscous tangential velocity C-------------------------------------------------------------- INCLUDE 'XFOIL.INC' C DO 10 IS=1, 2 UEDG(1,IS) = 0. DO 110 IBL=2, NBL(IS) I = IPAN(IBL,IS) UEDG(IBL,IS) = VTI(IBL,IS)*QVIS(I) 110 CONTINUE 10 CONTINUE C RETURN END SUBROUTINE QVFUE C-------------------------------------------------------------- C Sets panel viscous tangential velocity from viscous Ue C-------------------------------------------------------------- INCLUDE 'XFOIL.INC' C DO 1 IS=1, 2 DO 10 IBL=2, NBL(IS) I = IPAN(IBL,IS) QVIS(I) = VTI(IBL,IS)*UEDG(IBL,IS) 10 CONTINUE 1 CONTINUE C RETURN END SUBROUTINE QISET C------------------------------------------------------- C Sets inviscid panel tangential velocity for C current alpha. C------------------------------------------------------- INCLUDE 'XFOIL.INC' C COSA = COS(ALFA) SINA = SIN(ALFA) C DO 5 I=1, N+NW QINV (I) = COSA*QINVU(I,1) + SINA*QINVU(I,2) QINV_A(I) = -SINA*QINVU(I,1) + COSA*QINVU(I,2) 5 CONTINUE C RETURN END SUBROUTINE GAMQV INCLUDE 'XFOIL.INC' C DO 10 I=1, N GAM(I) = QVIS(I) GAM_A(I) = QINV_A(I) 10 CONTINUE C RETURN END SUBROUTINE STMOVE C--------------------------------------------------- C Moves stagnation point location to new panel. C--------------------------------------------------- INCLUDE 'XFOIL.INC' C C---- locate new stagnation point arc length SST from GAM distribution ISTOLD = IST CALL STFIND C IF(ISTOLD.EQ.IST) THEN C C----- recalculate new arc length array CALL XICALC C ELSE C CCC WRITE(*,*) 'STMOVE: Resetting stagnation point' C C----- set new BL position -> panel position pointers CALL IBLPAN C C----- set new inviscid BL edge velocity UINV from QINV CALL UICALC C C----- recalculate new arc length array CALL XICALC C C----- set BL position -> system line pointers CALL IBLSYS C IF(IST.GT.ISTOLD) THEN C------ increase in number of points on top side (IS=1) IDIF = IST-ISTOLD C ITRAN(1) = ITRAN(1) + IDIF ITRAN(2) = ITRAN(2) - IDIF C C------ move top side BL variables downstream DO 110 IBL=NBL(1), IDIF+2, -1 CTAU(IBL,1) = CTAU(IBL-IDIF,1) THET(IBL,1) = THET(IBL-IDIF,1) DSTR(IBL,1) = DSTR(IBL-IDIF,1) UEDG(IBL,1) = UEDG(IBL-IDIF,1) 110 CONTINUE C C------ set BL variables between old and new stagnation point DUDX = UEDG(IDIF+2,1)/XSSI(IDIF+2,1) DO 115 IBL=IDIF+1, 2, -1 CTAU(IBL,1) = CTAU(IDIF+2,1) THET(IBL,1) = THET(IDIF+2,1) DSTR(IBL,1) = DSTR(IDIF+2,1) UEDG(IBL,1) = DUDX * XSSI(IBL,1) 115 CONTINUE C C------ move bottom side BL variables upstream DO 120 IBL=2, NBL(2) CTAU(IBL,2) = CTAU(IBL+IDIF,2) THET(IBL,2) = THET(IBL+IDIF,2) DSTR(IBL,2) = DSTR(IBL+IDIF,2) UEDG(IBL,2) = UEDG(IBL+IDIF,2) 120 CONTINUE C ELSE C------ increase in number of points on bottom side (IS=2) IDIF = ISTOLD-IST C ITRAN(1) = ITRAN(1) - IDIF ITRAN(2) = ITRAN(2) + IDIF C C------ move bottom side BL variables downstream DO 210 IBL=NBL(2), IDIF+2, -1 CTAU(IBL,2) = CTAU(IBL-IDIF,2) THET(IBL,2) = THET(IBL-IDIF,2) DSTR(IBL,2) = DSTR(IBL-IDIF,2) UEDG(IBL,2) = UEDG(IBL-IDIF,2) 210 CONTINUE C C------ set BL variables between old and new stagnation point DUDX = UEDG(IDIF+2,2)/XSSI(IDIF+2,2) c write(*,*) 'idif Ue xi dudx', c & idif, UEDG(idif+2,2), xssi(idif+2,2), dudx DO 215 IBL=IDIF+1, 2, -1 CTAU(IBL,2) = CTAU(IDIF+2,2) THET(IBL,2) = THET(IDIF+2,2) DSTR(IBL,2) = DSTR(IDIF+2,2) UEDG(IBL,2) = DUDX * XSSI(IBL,2) 215 CONTINUE c write(*,*) 'Uenew xinew', idif+1, uedg(idif+1,2), xssi(idif+1,2) C C------ move top side BL variables upstream DO 220 IBL=2, NBL(1) CTAU(IBL,1) = CTAU(IBL+IDIF,1) THET(IBL,1) = THET(IBL+IDIF,1) DSTR(IBL,1) = DSTR(IBL+IDIF,1) UEDG(IBL,1) = UEDG(IBL+IDIF,1) 220 CONTINUE ENDIF C C----- tweak Ue so it's not zero, in case stag. point is right on node UEPS = 1.0E-7 DO IS = 1, 2 DO IBL = 2, NBL(IS) I = IPAN(IBL,IS) IF(UEDG(IBL,IS).LE.UEPS) THEN UEDG(IBL,IS) = UEPS QVIS(I) = VTI(IBL,IS)*UEPS GAM(I) = VTI(IBL,IS)*UEPS ENDIF ENDDO ENDDO C ENDIF C C---- set new mass array since Ue has been tweaked DO 50 IS=1, 2 DO 510 IBL=2, NBL(IS) MASS(IBL,IS) = DSTR(IBL,IS)*UEDG(IBL,IS) 510 CONTINUE 50 CONTINUE C RETURN END SUBROUTINE UESET C--------------------------------------------------------- C Sets Ue from inviscid Ue plus all source influence C--------------------------------------------------------- INCLUDE 'XFOIL.INC' C DO 1 IS=1, 2 DO 10 IBL=2, NBL(IS) I = IPAN(IBL,IS) C DUI = 0. DO 100 JS=1, 2 DO 1000 JBL=2, NBL(JS) J = IPAN(JBL,JS) UE_M = -VTI(IBL,IS)*VTI(JBL,JS)*DIJ(I,J) DUI = DUI + UE_M*MASS(JBL,JS) 1000 CONTINUE 100 CONTINUE C UEDG(IBL,IS) = UINV(IBL,IS) + DUI C 10 CONTINUE 1 CONTINUE C RETURN END SUBROUTINE DSSET INCLUDE 'XFOIL.INC' C DO 1 IS=1, 2 DO 10 IBL=2, NBL(IS) DSTR(IBL,IS) = MASS(IBL,IS) / UEDG(IBL,IS) 10 CONTINUE 1 CONTINUE C RETURN END