NACA0012라는 에어포일을 PANEL METHOD라는 방법으로 분석하는 코딩을 해봤습니다.

포트란이 아닌 전산유체시뮬레이션으로 돌렸을 때는 CL값이 어느정도 증가하다가 곡선을 이루며 감소하는데

이 코드로써 얻은 각도(AOA)에 따른 양력계수(CL)값의 그래프가 직선으로 나오네요...직선이 되면 안되는데
(input=Angle of Attack/ output=CL)

어느부분을 잘못한 걸까요..... 

 

 

 

 

C ==============================================================
C VORTEX PANEL METHOD
C ==============================================================
      PROGRAM MAIN
      COMMON/GEO/XE(200),YE(200),XC(200),YC(200),DL(200),XT(200)
     &        ,YT(200)
      COMMON/DATA/N,NP1,ALPA,PI,PI2,GAMMA
      COMMON/VEL/UIJ,VIJ
      COMMON/SIM/A(200,200)
      COMMON/FORCE/GAM(200),CP(200)
      DIMENSION RHS(200)
C --------------------------------------------
      OPEN(6,FILE='OUT.OUT',STATUS='unknown')
      PI=ATAN(1.)*4.
      PI2=PI*2.
C --------------------------------------------------------------
C INPUT Angle of Attack
C --------------------------------------------------------------

      WRITE(*,*)' ALPHA(deg.)=? '
      READ(*,*)ALPA
      ALPA=ALPA*PI/180.
C --------------------------------------------------------------
C SET THE AIRFOIL SHAPE
      CALL PROFILE
C --------------------------------------------------------------

C CALCULATE CONTROL POINTS OF ELEMENTS and TANGENTIAL DIRECTION
      DO 30 I=1,N
      IP=I+1
      IF(I.EQ.N)IP=1
      XC(I)=(XE(I)+XE(IP))/2.
      YC(I)=(YE(I)+YE(IP))/2.
      DX=XE(IP)-XE(I)
      DY=YE(IP)-YE(I)
      DL(I)=SQRT(DX**2+DY**2)  ! LENGTH OF ELEMENTS
      XT(I)=DX/DL(I)   ! COS(PHI_I)
      YT(I)=DY/DL(I)   ! SIN(PHI_I)
30    CONTINUE
C --------------------------------------------------------------
C CALCULATE OEEFFICIENT MATRIX
      DO 60 I=1,N
      DO 60 J=1,N
      JP=J+1
      IF(I.NE.J)THEN
      CALL VTCVEL(XC(I),YC(I),XE(J),YE(J),XE(JP),YE(JP),XT(J),YT(J))
      A(I,J)=(XT(I)*UIJ+YT(I)*VIJ)/PI2
      ELSE
      A(I,J)=-0.5
      ENDIF
60    CONTINUE
C --------------------------------------------------------------
C APPLY THE KUTTA CONDITION
      DO J=1,N
      A(N,J)=0.0
      ENDDO
      A(N,2)=1.0
      A(N,N-1)=1.0
C --------------------------------------------------------------
C INVERSE COEFFICIENT MATRIX
C --------------------------------------------------------------
      CALL SIMUL(N)
      IF(N.EQ.0)STOP
C --------------------------------------------------------------
C CALCULATER RIGHT HAND VECTOR
      UINF=COS(ALPA)
      VINF=SIN(ALPA)
      DO 90 I=1,N
      RHS(I)=-(UINF*XT(I)+VINF*YT(I))
90    CONTINUE
      RHS(N)=0.0
C --------------------------------------------------------------
      DO 110 I=1,N
      GAM(I)=0.0
      DO 110 J=1,N
      GAM(I)=GAM(I)+A(I,J)*RHS(J)
110   CONTINUE
C --------------------------------------------------------------
C CALCULATE AERODYNAMIC VALUES
C --------------------------------------------------------------
      CM=0.0
      DO 80 I=1,N
      CP(I)=1.-GAM(I)**2
      CM=CM-CP(I)*DL(I)*(XT(I)*(XC(I)-0.25)+YT(I)*YC(I))
80    CONTINUE
C --------------------------------------------------------------
C CALCULATE VORTEX STRENGTH
      GAMMA=0.0
      DO 150 I=1,N
150    GAMMA=GAMMA+DL(I)*GAM(I)
      CL=-GAMMA*2.4
      WRITE(*,560)CL,CM
C --------------------------------------------------------------
C WRITE OUT THE RESULTS
      WRITE(6,500)
500   FORMAT(2X,70('-')/5X,'XE(I)',5X,'YE(I)',5X,'XC(I)',5X,'YC(I)'
     &        ,4X,'RHS(I)',6X,'V(I)',5X,'CP(I)'/2X,70('-'))
      DO I=1,NP1
      ENDDO
      WRITE(6,510)XE(I),YE(I),XC(I),YC(I),RHS(I),GAM(I),CP(I)
510   FORMAT(2X,7F10.4)
      WRITE(6,560)CL,CM
560   FORMAT(//5X,'LIFT COEEF.=',F8.4,'PTCH MOMENT COEEF.=',F8.4)
C --------------------------------------------------------------
      PAUSE 
      END
C ===================================================================================================
      SUBROUTINE PROFILE
C ---------------------------------------------------------------------------------------------------
C SUBPROGRAM FOR THE AIRFOIL SURFACE COORDINATES GENERATION.
      COMMON/GEO/XE(200),YE(200),XC(200),YC(200),DL(200),XT(200)
     &         ,YT(200)
      COMMON/DATA/N,NP1,ALPA,PI,PI2,GAMMA
      DIMENSION YCAMB(100),XCAMB(100)
C -----------------------------------------
C THICKNESS DISTRIBUTION OF NACA-0012 AIRFOIL
      THICK(X,TMAX)=TMAX/0.2*(0.296*SQRT(X)-0.126*X-0.3516*X**2
     &          +0.2843*X**3-0.1015*X**4)
C Where the "TMAX" is the max. thickness of the NACA 4 or 5 digit airfoil
C -----------------------------------------------------------------------------------------------------
C FOR NACA-0012 AIRFOIL
   CMB=0.0     ! MAX. CAMBER
      PCM=0     ! POSITION OF MAX. CAMBER
      TMAX=0.12     ! MAX. THICKNESS
C -----------------------------------------------------------------------------------------------------
CAMBER COORDINATES OF AIRFOIL
      M=60
      N=2*M
      NP1=N+1
      DELTHE=0.5*PI/FLOAT(M)
      XCAMB(1)=0.0
      YCAMB(1)=0.0
      THE=PI
      DO I=2,M
      THE=THE-DELTHE
      XX=(1.+COS(THE))
      XCAMB(I)=XX
      IF(XX.LE.PCM)THEN
      YCAMB(I)=CMB/PCM**2*(2.*PCM*XX-XX**2)
      ELSE IF(XX.GT.PCM) THEN
      YCAMB(I)=CMB/(1.-PCM)**2*(1.0+2.*PCM*(XX-1.)-XX**2)
      ENDIF
      ENDDO
      XCAMB(M+1)=1.0
      YCAMB(M+1)=0.0
C -----------------------------------------------------------------------------------------------------
C AIRFOIL SURFACE COORDINATES
      K=1
      L=NP1
      XE(K)=1.0
      YE(K)=0.0
      XE(L)=1.0
      YE(L)=0.0
      DO I=M,2,-1
      K=K+1
      L=L-1
      XX=XCAMB(I)
      IP1=I+1
      IM1=I-1
      DX=XCAMB(IP1)-XCAMB(IM1)
      DY=YCAMB(IP1)-YCAMB(IM1)
      DLL=SQRT(DX**2+DY**2)
      SN=DY/DLL
      CS=DX/DLL
      DELTA=THICK(XX,TMAX)
      XE(K)=XCAMB(I)-DELTA*SN
      YE(K)=YCAMB(I)+DELTA*CS
      XE(L)=XCAMB(I)+DELTA*SN
      YE(L)=YCAMB(I)-DELTA*CS
      ENDDO
      XE(M+1)=0.0
      YE(M+1)=0.0
      RETURN
      END
C ===================================================================================================
C SUBPROGRAM FOR LINEAR EQUATION
C -----------------------------------------------------------------------------------------------------
      SUBROUTINE SIMUL(N)
      COMMON/SIM/A(200,200)
      DO 18 K=1,N
      PIVOT=A(K,K)
      IF(ABS(PIVOT).GT.0.00000000001) GOTO 13
      GOTO 30
13    DO 14 J=1,N
14    A(K,J)=A(K,J)/PIVOT
   A(K,K)=1./PIVOT
      DO 18 I=1,N
      AIJCK=A(I,K)
      IF(I.EQ.K) GOTO 18
      A(I,K)=-AIJCK/PIVOT
      DO 17 J=1,N
17    IF(J.NE.K) A(I,J)=A(I,J)-AIJCK*A(K,J)
18    CONTINUE
   RETURN
30    N=0
   RETURN
      END
C ===================================================================================================
      SUBROUTINE VTCVEL(XO,YO,X1,Y1,X2,Y2,XT,YT)
C ===================================================================================================
C SUBPROGRAM FOR VELOCITY VECTOR
    COMMON/VEL/UX,VY
C -------------------------------------
C COORDINATE TRANSFORM
      XN1=XO-X1
      XN2=XO-X2
      YN1=YO-Y1
      YN2=YO-Y2
      RC1=XN1*XT+YN1*YT    ! X'_i
      RC2=XN2*XT+YN2*YT    ! X'_i+1
      RS1=YN1*XT-XN1*YT    ! Y'_i
      R1=XN1**2+YN1**2    ! r_i,j
      R2=XN2**2+YN2**2    ! r_i,j+1
      TH1=ATAN2(RS1,RC1)    ! THETA_i,j
      TH2=ATAN2(RS1,RC2)    ! THETA_i,j+1
      TH2N1=TH2-TH1
      ALR1O2=0.5*ALOG(R2/R1)
      U1=-TH2N1
      V1=-ALR1O2
C INVERSE TRANSFORM FOR VELOCITY VEXTOR
      UX=U1*XT-V1*YT
      VY=U1*YT+V1*XT
      RETURN
      END
C ===================================================================================================