<div>NACA0012라는 에어포일을 PANEL METHOD라는 방법으로 분석하는 코딩을 해봤습니다.</div> <div>포트란이 아닌 전산유체시뮬레이션으로 돌렸을 때는 CL값이 어느정도 증가하다가 곡선을 이루며 감소하는데</div> <div>이 코드로써 얻은 각도(AOA)에 따른 양력계수(CL)값의 그래프가 직선으로 나오네요...직선이 되면 안되는데 <br>(input=Angle of Attack/ output=CL)</div> <div>어느부분을 잘못한 걸까요..... </div> <div> </div> <div> </div> <div> </div> <div> </div> <div>C ==============================================================<br>C VORTEX PANEL METHOD<br>C ==============================================================<br>      PROGRAM MAIN<br>      COMMON/GEO/XE(200),YE(200),XC(200),YC(200),DL(200),XT(200)<br>     &        ,YT(200)<br>      COMMON/DATA/N,NP1,ALPA,PI,PI2,GAMMA<br>      COMMON/VEL/UIJ,VIJ<br>      COMMON/SIM/A(200,200)<br>      COMMON/FORCE/GAM(200),CP(200)<br>      DIMENSION RHS(200)<br>C --------------------------------------------<br>      OPEN(6,FILE='OUT.OUT',STATUS='unknown')<br>      PI=ATAN(1.)*4.<br>      PI2=PI*2.<br>C --------------------------------------------------------------<br>C INPUT Angle of Attack<br>C --------------------------------------------------------------</div> <div>      WRITE(*,*)' ALPHA(deg.)=? '<br>      READ(*,*)ALPA<br>      ALPA=ALPA*PI/180.<br>C --------------------------------------------------------------<br>C SET THE AIRFOIL SHAPE<br>      CALL PROFILE<br>C --------------------------------------------------------------</div> <div>C CALCULATE CONTROL POINTS OF ELEMENTS and TANGENTIAL DIRECTION<br>      DO 30 I=1,N<br>      IP=I+1<br>      IF(I.EQ.N)IP=1<br>      XC(I)=(XE(I)+XE(IP))/2.<br>      YC(I)=(YE(I)+YE(IP))/2.<br>      DX=XE(IP)-XE(I)<br>      DY=YE(IP)-YE(I)<br>      DL(I)=SQRT(DX**2+DY**2)  ! LENGTH OF ELEMENTS<br>      XT(I)=DX/DL(I)   ! COS(PHI_I)<br>      YT(I)=DY/DL(I)   ! SIN(PHI_I)<br>30    CONTINUE<br>C --------------------------------------------------------------<br>C CALCULATE OEEFFICIENT MATRIX<br>      DO 60 I=1,N<br>      DO 60 J=1,N<br>      JP=J+1<br>      IF(I.NE.J)THEN<br>      CALL VTCVEL(XC(I),YC(I),XE(J),YE(J),XE(JP),YE(JP),XT(J),YT(J))<br>      A(I,J)=(XT(I)*UIJ+YT(I)*VIJ)/PI2<br>      ELSE<br>      A(I,J)=-0.5<br>      ENDIF<br>60    CONTINUE<br>C --------------------------------------------------------------<br>C APPLY THE KUTTA CONDITION<br>      DO J=1,N<br>      A(N,J)=0.0<br>      ENDDO<br>      A(N,2)=1.0<br>      A(N,N-1)=1.0<br>C --------------------------------------------------------------<br>C INVERSE COEFFICIENT MATRIX<br>C --------------------------------------------------------------<br>      CALL SIMUL(N)<br>      IF(N.EQ.0)STOP<br>C --------------------------------------------------------------<br>C CALCULATER RIGHT HAND VECTOR<br>      UINF=COS(ALPA)<br>      VINF=SIN(ALPA)<br>      DO 90 I=1,N<br>      RHS(I)=-(UINF*XT(I)+VINF*YT(I))<br>90    CONTINUE<br>      RHS(N)=0.0<br>C --------------------------------------------------------------<br>      DO 110 I=1,N<br>      GAM(I)=0.0<br>      DO 110 J=1,N<br>      GAM(I)=GAM(I)+A(I,J)*RHS(J)<br>110   CONTINUE<br>C --------------------------------------------------------------<br>C CALCULATE AERODYNAMIC VALUES<br>C --------------------------------------------------------------<br>      CM=0.0<br>      DO 80 I=1,N<br>      CP(I)=1.-GAM(I)**2<br>      CM=CM-CP(I)*DL(I)*(XT(I)*(XC(I)-0.25)+YT(I)*YC(I))<br>80    CONTINUE<br>C --------------------------------------------------------------<br>C CALCULATE VORTEX STRENGTH<br>      GAMMA=0.0<br>      DO 150 I=1,N<br>150    GAMMA=GAMMA+DL(I)*GAM(I)<br>      CL=-GAMMA*2.4<br>      WRITE(*,560)CL,CM<br>C --------------------------------------------------------------<br>C WRITE OUT THE RESULTS<br>      WRITE(6,500)<br>500   FORMAT(2X,70('-')/5X,'XE(I)',5X,'YE(I)',5X,'XC(I)',5X,'YC(I)'<br>     &        ,4X,'RHS(I)',6X,'V(I)',5X,'CP(I)'/2X,70('-'))<br>      DO I=1,NP1<br>      ENDDO<br>      WRITE(6,510)XE(I),YE(I),XC(I),YC(I),RHS(I),GAM(I),CP(I)<br>510   FORMAT(2X,7F10.4)<br>      WRITE(6,560)CL,CM<br>560   FORMAT(//5X,'LIFT COEEF.=',F8.4,'PTCH MOMENT COEEF.=',F8.4)<br>C --------------------------------------------------------------<br>      PAUSE <br>      END<br>C ===================================================================================================<br>      SUBROUTINE PROFILE<br>C ---------------------------------------------------------------------------------------------------<br>C SUBPROGRAM FOR THE AIRFOIL SURFACE COORDINATES GENERATION.<br>      COMMON/GEO/XE(200),YE(200),XC(200),YC(200),DL(200),XT(200)<br>     &         ,YT(200)<br>      COMMON/DATA/N,NP1,ALPA,PI,PI2,GAMMA<br>      DIMENSION YCAMB(100),XCAMB(100)<br>C -----------------------------------------<br>C THICKNESS DISTRIBUTION OF NACA-0012 AIRFOIL<br>      THICK(X,TMAX)=TMAX/0.2*(0.296*SQRT(X)-0.126*X-0.3516*X**2<br>     &          +0.2843*X**3-0.1015*X**4)<br>C Where the "TMAX" is the max. thickness of the NACA 4 or 5 digit airfoil<br>C -----------------------------------------------------------------------------------------------------<br>C FOR NACA-0012 AIRFOIL<br>   CMB=0.0     ! MAX. CAMBER<br>      PCM=0     ! POSITION OF MAX. CAMBER<br>      TMAX=0.12     ! MAX. THICKNESS<br>C -----------------------------------------------------------------------------------------------------<br>CAMBER COORDINATES OF AIRFOIL<br>      M=60<br>      N=2*M<br>      NP1=N+1<br>      DELTHE=0.5*PI/FLOAT(M)<br>      XCAMB(1)=0.0<br>      YCAMB(1)=0.0<br>      THE=PI<br>      DO I=2,M<br>      THE=THE-DELTHE<br>      XX=(1.+COS(THE))<br>      XCAMB(I)=XX<br>      IF(XX.LE.PCM)THEN<br>      YCAMB(I)=CMB/PCM**2*(2.*PCM*XX-XX**2)<br>      ELSE IF(XX.GT.PCM) THEN<br>      YCAMB(I)=CMB/(1.-PCM)**2*(1.0+2.*PCM*(XX-1.)-XX**2)<br>      ENDIF<br>      ENDDO<br>      XCAMB(M+1)=1.0<br>      YCAMB(M+1)=0.0<br>C -----------------------------------------------------------------------------------------------------<br>C AIRFOIL SURFACE COORDINATES<br>      K=1<br>      L=NP1<br>      XE(K)=1.0<br>      YE(K)=0.0<br>      XE(L)=1.0<br>      YE(L)=0.0<br>      DO I=M,2,-1<br>      K=K+1<br>      L=L-1<br>      XX=XCAMB(I)<br>      IP1=I+1<br>      IM1=I-1<br>      DX=XCAMB(IP1)-XCAMB(IM1)<br>      DY=YCAMB(IP1)-YCAMB(IM1)<br>      DLL=SQRT(DX**2+DY**2)<br>      SN=DY/DLL<br>      CS=DX/DLL<br>      DELTA=THICK(XX,TMAX)<br>      XE(K)=XCAMB(I)-DELTA*SN<br>      YE(K)=YCAMB(I)+DELTA*CS<br>      XE(L)=XCAMB(I)+DELTA*SN<br>      YE(L)=YCAMB(I)-DELTA*CS<br>      ENDDO<br>      XE(M+1)=0.0<br>      YE(M+1)=0.0<br>      RETURN<br>      END<br>C ===================================================================================================<br>C SUBPROGRAM FOR LINEAR EQUATION<br>C -----------------------------------------------------------------------------------------------------<br>      SUBROUTINE SIMUL(N)<br>      COMMON/SIM/A(200,200)<br>      DO 18 K=1,N<br>      PIVOT=A(K,K)<br>      IF(ABS(PIVOT).GT.0.00000000001) GOTO 13<br>      GOTO 30<br>13    DO 14 J=1,N<br>14    A(K,J)=A(K,J)/PIVOT<br>   A(K,K)=1./PIVOT<br>      DO 18 I=1,N<br>      AIJCK=A(I,K)<br>      IF(I.EQ.K) GOTO 18<br>      A(I,K)=-AIJCK/PIVOT<br>      DO 17 J=1,N<br>17    IF(J.NE.K) A(I,J)=A(I,J)-AIJCK*A(K,J)<br>18    CONTINUE<br>   RETURN<br>30    N=0<br>   RETURN<br>      END<br>C ===================================================================================================<br>      SUBROUTINE VTCVEL(XO,YO,X1,Y1,X2,Y2,XT,YT)<br>C ===================================================================================================<br>C SUBPROGRAM FOR VELOCITY VECTOR<br>    COMMON/VEL/UX,VY<br>C -------------------------------------<br>C COORDINATE TRANSFORM<br>      XN1=XO-X1<br>      XN2=XO-X2<br>      YN1=YO-Y1<br>      YN2=YO-Y2<br>      RC1=XN1*XT+YN1*YT    ! X'_i<br>      RC2=XN2*XT+YN2*YT    ! X'_i+1<br>      RS1=YN1*XT-XN1*YT    ! Y'_i<br>      R1=XN1**2+YN1**2    ! r_i,j<br>      R2=XN2**2+YN2**2    ! r_i,j+1<br>      TH1=ATAN2(RS1,RC1)    ! THETA_i,j<br>      TH2=ATAN2(RS1,RC2)    ! THETA_i,j+1<br>      TH2N1=TH2-TH1<br>      ALR1O2=0.5*ALOG(R2/R1)<br>      U1=-TH2N1<br>      V1=-ALR1O2<br>C INVERSE TRANSFORM FOR VELOCITY VEXTOR<br>      UX=U1*XT-V1*YT<br>      VY=U1*YT+V1*XT<br>      RETURN<br>      END<br>C ===================================================================================================</div>