PROGRAM VLM
use MSIMSL
IMPLICIT NONE
REAL*8,allocatable :: XB(:),A(:,:),BL(:),SOLV(:),PHI(:),PHIW(:),PHIOLD(:)
REAL*8,allocatable :: CF(:)
INTEGER :: I,J,NW,NB,IT
REAL*8 :: UINF,ALPHA,PI,DT,CL,AINF,DX,TIME,T,T0,MACH,SUM1
OPEN (6,FILE="CLT.DAT")
OPEN (4,FILE="CIRCUL.DAT")
PI = Dacos(-1.0D0)
ALPHA = 2.0D0 * PI / 180.0d0
DX = 0.01D0
MACH = 0.000001D0
! AINF = 300.0d0
! UINF = MACH * AINF
UINF = 1.0D0
AINF = UINF / MACH
DT = 0.01d0 !1.5D0 / ( 1.0d0/MACH+1.0d0)
TIME = 10.0D0
NW = TIME / DT
NB = 1.0D0 / DX + 1
ALLOCATE(XB(NB),A(NB,NB),BL(NB),SOLV(NB),PHI(NB),PHIOLD(NB),PHIW(NW),CF(NB-1))
CALL READMESH (NB,DX,XB)
PHI = 0.0D0 ; PHIW = 0.0D0 ; PHIOLD =0.0D0
SUM1 = 0.0D0
DO IT = 1 , TIME/DT
T= IT * DT
T0 = 0.0D0
PHIW(IT) = PHI(NB)
PHIOLD = PHI
CALL COEFF_MATRIX (DT,T,T0,UINF,AINF,NB,XB,A)
CALL RHS_VECTOR (NB,NW,IT,SUM1,XB,PHIW,ALPHA,UINF,AINF,DT,T,T0,BL)
CALL DLSARG (NB, A, NB, BL, 1, SOLV)
DO I = 1 , NB
PHI(I) = SOLV(I)
END DO
CALL PRESS(DX,DT,UINF,NB,PHI,PHIOLD,CF,CL)
WRITE(6,*)IT*DT , CL/(2.0d0*PI*alpha)
WRITE(*,*)IT*DT , CL/(2.0d0*PI*alpha)
END DO ! TIME
DO I = NB-1 , 1 ,-1
WRITE(4,102)XB(I) , PHI(I)
END DO
102 Format (1x,10(1x,f14.5))
STOP
END
!-----------------------------------------------------------------------------------------
!-------------------------------- SUBROUTINES --------------------------------------------
!-----------------------------------------------------------------------------------------
SUBROUTINE PRESS(DX,DT,UINF,NB,PHI,PHIOLD,CF,CL)
Implicit NONE
INTEGER,INTENT(IN) :: NB
REAL*8,INTENT(IN) :: PHI(NB),PHIOLD(NB),DT,UINF,DX
REAL*8,INTENT(OUT) :: CF(NB-1),CL
REAL*8 :: SUM
INTEGER :: I,J
CF = 0.0D0
DO I = 1 , NB-1
SUM = 0.0D0
DO J = 1 , I
SUM = SUM + (PHI(J) - PHIOLD(J) )
END DO
CF(I) = 1.0D0 * UINF * PHI(I) + 1.0D0 / DT * SUM * DX
END DO
CL = 0.0D0
DO I = 1 , NB-1
CL = CL + CF(i) / (0.5 * 1.0D0 * UINF**2 * 1.0D0)
END DO
END SUBROUTINE PRESS
SUBROUTINE RHS_VECTOR (NB,NW,IT,SUM1,XB,PHIW,ALPHA,UINF,AINF,DT,T,T0,BL)
Implicit NONE
INTEGER,INTENT(IN) :: NB,NW,IT
REAL*8,INTENT(IN) :: XB(NB),PHIW(NW),ALPHA,UINF,AINF,DT,T,T0
REAL*8,INTENT(INOUT) :: SUM1
REAL*8,INTENT(OUT) :: BL(NB)
REAL*8 :: X,PI,PI2,DX,DTL,SUM2
INTEGER::I,J,K
PI = Dacos(-1.0D0)
PI2 = 2.0D0*PI
SUM1 = SUM1 + PHIW(IT)
BL(NB) = -SUM1
Do I = 1 , NB-1 ! loop over collocation points , X0 : Collocation Point
X = ( 3.0D0*XB(I) + 1.0D0*XB(I+1) ) / 4.0D0
SUM2 = 0.0D0
DO J = 1 , IT
DTL = (IT+2-J) * DT
DX = X - (1.0D0 + UINF * DTL)
! if ((AINF * DTL)**2 - DX**2 .lt. 0.0d0) then
! write(*,*)'pause3'
! write(*,*)t ,i ,j , x , 1.0D0 + UINF * DT
! end if
SUM2 = SUM2 - PHIW(J) / (PI2 * AINF * DTL * DX) * DSQRT ( (AINF * DTL)**2 - DX**2)
! SUM2 = SUM2 - PHIW(J) / (PI2 * DX)
END DO
BL(I) = -UINF * alpha - SUM2
END DO
END SUBROUTINE RHS_VECTOR
SUBROUTINE COEFF_MATRIX (DT,T,T0,UINF,AINF,NB,XB,A)
Implicit NONE
INTEGER,INTENT(IN) :: NB
REAL*8,INTENT(IN) :: XB(NB),DT,T,T0,UINF,AINF
REAL*8,INTENT(OUT) :: A(NB,NB)
REAL*8 :: X,X0,G,PI,PI2,DX,DTL
INTEGER::I,J,K
PI = Dacos(-1.0D0)
PI2 = 2.0D0*PI
DTL = T - T
A=0.0D0
Do I=1,NB-1 ! loop over collocation points , X : Collocation Point
X = ( 3.0D0*XB(I) + 1.0D0*XB(I+1) ) / 4.0D0
Do J=1,NB-1 ! loop over ELEMENTS OF AIRFOIL (B)
X0 = ( 1.0D0*XB(J) + 3.0D0*XB(J+1)) / 4.0D0
DX = X - X0
! if ((AINF*DT)**2 - (DX-UINF*DT)**2 .lt. 0.0d0) then
! write(*,*)'pause1'
! write(*,*) i , j , x , DX , (AINF*DT)**2 - (DX-UINF*DT)**2
! end if
G = -1.0D0 / (PI2 * AINF *DX * DT) * DSQRT ((AINF*DT)**2 - (DX-UINF*DT)**2)
A(I,J) = G
END DO ! J
DX = X - ( 1.0D0 + UINF * DT )
! if ((AINF * DT)**2 - DX**2 .lt. 0.0d0) then
! write(*,*)'pause2'
! write(*,*)t ,i , x , (AINF * DT)**2 , DX**2
! end if
G = -1.0D0/(PI2 * AINF * DT * DX) * DSQRT ((AINF*DT)**2 - DX **2)
! G = -1.0D0/(PI2 * DX)
A(I,NB) = G
END DO ! I
A(NB,:) = 1.0D0
END SUBROUTINE COEFF_MATRIX
SUBROUTINE READMESH (NB,DX,XB)
Implicit NONE
INTEGER,INTENT(IN):: NB
REAL*8 , INTENT (IN):: DX
REAL*8 , INTENT (OUT):: XB(NB)
INTEGER::I,J
DO I = 1 , NB
XB(I) = 1.0D0 - (I-1)*DX
END DO
END SUBROUTINE READMESH