[sezo85]

Hi everyone,

I am a fresher in numerical techniques and I am developing a code to simulate a compressor. I am using FORTRAN 95, my code will be a quasi-1D solver. I am using common Euler equations, with an additional equation to solve the tangential momentum.
Till now I have been using the MacCormark scheme for the time integration together with a finite difference approach.
Now I am stuck since I am dealing with a compressor characterized by a very long domain and very low mach numbers...to fix the convergence problem I implemented Crank Nicolson ..and for the moment it doesn't want to work..
Could anyone give me any advice??? To implement this scheme I followed the book by Anderson..I reported my code below..

ADDITIONAL INFORMATION:
-my equations are written in the conservative form;
-DF are the derivatives of the fluxes respect to U;
-1,2,3,4 correspond to mass conservation, axial momentum, tangential momentum and energy equation;
-S terms collect all the sources to simulate a compressor and changes in cross section;
-D is what is know at every timestep..as explained in Anderson book.

Please If anyone can spot an error..it would be great

Thank you very much.

SUBROUTINE CrankNicolson

USE Var_Declar
USE BASIS_R
IMPLICIT NONE

DT(PP)=0.0000005d0

!ARTIFICIAL VISCOSITY
DO j=2,SIZE_ST-1

VIS1(j,PP)=0.8D0*ABS(P(j+1,PP)-2.0D0*P(j,PP)+P(j-1,PP))/(P(j+1,PP)+2.0D0*P(j,PP)+P(j-1,PP))*(U1(j+1,PP)-2.0D0*U1(j,PP)+U1(j-1,PP))
VIS2(j,PP)=0.8D0*ABS(P(j+1,PP)-2.0D0*P(j,PP)+P(j-1,PP))/(P(j+1,PP)+2.0D0*P(j,PP)+P(j-1,PP))*(U2(j+1,PP)-2.0D0*U2(j,PP)+U2(j-1,PP))
VIS3(j,PP)=0.8D0*ABS(P(j+1,PP)-2.0D0*P(j,PP)+P(j-1,PP))/(P(j+1,PP)+2.0D0*P(j,PP)+P(j-1,PP))*(U3(j+1,PP)-2.0D0*U3(j,PP)+U3(j-1,PP))
VIS4(j,PP)=0.8D0*ABS(P(j+1,PP)-2.0D0*P(j,PP)+P(j-1,PP))/(P(j+1,PP)+2.0D0*P(j,PP)+P(j-1,PP))*(U4(j+1,PP)-2.0D0*U4(j,PP)+U4(j-1,PP))

END DO

DO j=1,SIZE_ST
!Evaluation DF matrix of known terms at the previous timestep
!DF_1(j,PP)=(F1(J+1,pp)-F1(J-1,pp))/(2*(U1(J+1,pp)-U1(J-1,pp)))!0.0D0
!
!DF_2(j,PP)=(F2(J+1,pp)-F2(J-1,pp))/(2*(U2(J+1,pp)-U2(J-1,pp)))!(3-GAMMA)*U2(j,PP)/U1(j,PP)
!
!DF_3(j,PP)=(F3(J+1,pp)-F3(J-1,pp))/(2*(U3(J+1,pp)-U3(J-1,pp)))!(U2(j,PP)/U1(j,PP))
!
!DF_4(j,PP)=(F4(J+1,pp)-F4(J-1,pp))/(2*(U4(J+1,pp)-U4(J-1,pp)))!U2(j,PP)/U1(j,PP)*GAMMA

DF_1(j,PP)=0.0D0

DF_2(j,PP)=(3-GAMMA)*U2(j,PP)/U1(j,PP)

DF_3(j,PP)=(U2(j,PP)/U1(j,PP))

DF_4(j,PP)=U2(j,PP)/U1(j,PP)*GAMMA
ENDDO

DO j=2,SIZE_ST-1

!Evaluation source terms
S_1(j,PP)=0.0D0

S_2(j,PP)=J2_STAR(j,PP)*(LOG(SAREA(j+1,PP))-LOG(SAREA(j-1,PP)))/(2.0d0*DX(j,PP))+F_Z(j,PP)+J2_THETA(j,PP)*(MEAN_RA DIUS(j+1,PP)-MEAN_RADIUS(j-1,PP))/(2*DX(j,PP))

S_3(j,PP)=F_T(j,PP)/mean_radius(j,PP)+F_K(j,PP)*(MEAN_RADIUS(j+1,PP)-MEAN_RADIUS(j-1,PP))/(2.0d0*DX(j,PP))

S_4(j,PP)=W(j,PP)

! Evaluation vector D

D_1(j,PP)=U1(j,PP)-DT(PP)/(2.0d0*DX(j,PP))*(F1(J+1,PP)-F1(J-1,PP))+DT(PP)/(4.0d0*DX(j,PP))*(DF_1(j+1,PP)*U1(j+1,PP)-DF_1(j-1,PP)*U1(j-1,PP))-S_1(j,PP)*DT(PP)

D_2(j,PP)=U2(j,PP)-DT(PP)/(2.0d0*DX(j,PP))*(F2(J+1,PP)-F2(J-1,PP))+DT(PP)/(4.0d0*DX(j,PP))*(DF_2(j+1,PP)*U2(j+1,PP)-DF_2(j-1,PP)*U2(j-1,PP))-S_2(j,PP)*DT(PP)

D_3(j,PP)=U3(j,PP)-DT(PP)/(2.0d0*DX(j,PP))*(F3(J+1,PP)-F3(J-1,PP))+DT(PP)/(4.0d0*DX(j,PP))*(DF_3(j+1,PP)*U3(j+1,PP)-DF_3(j-1,PP)*U3(j-1,PP))-S_3(j,PP)*DT(PP)

D_4(j,PP)=U4(j,PP)-DT(PP)/(2.0d0*DX(j,PP))*(F4(J+1,PP)-F4(J-1,PP))+DT(PP)/(4.0d0*DX(j,PP))*(DF_4(j+1,PP)*U4(j+1,PP)-DF_4(j-1,PP)*U4(j-1,PP))-S_4(j,PP)*DT(PP)

END DO

DIAG_C(:,1)=DF_1(2:SIZE_ST-1,PP)*DT(PP)/(4.0d0*DX(1:SIZE_ST-2,PP))
DIAG_C(:,2)=DF_2(2:SIZE_ST-1,PP)*DT(PP)/(4.0d0*DX(1:SIZE_ST-2,PP))
DIAG_C(:,3)=DF_3(2:SIZE_ST-1,PP)*DT(PP)/(4.0d0*DX(1:SIZE_ST-2,PP))
DIAG_C(:,4)=DF_4(2:SIZE_ST-1,PP)*DT(PP)/(4.0d0*DX(1:SIZE_ST-2,PP))

DIAG_B(:,1)=1.0D0
DIAG_B(:,2)=1.0D0
DIAG_B(:,3)=1.0D0
DIAG_B(:,4)=1.0D0

DIAG_A(:,1)=-DF_1(2:SIZE_ST-1,PP)*DT(PP)/(4.0d0*DX(2:SIZE_ST-1,PP))
DIAG_A(:,2)=-DF_2(2:SIZE_ST-1,PP)*DT(PP)/(4.0d0*DX(2:SIZE_ST-1,PP))
DIAG_A(:,3)=-DF_3(2:SIZE_ST-1,PP)*DT(PP)/(4.0d0*DX(2:SIZE_ST-1,PP))
DIAG_A(:,4)=-DF_4(2:SIZE_ST-1,PP)*DT(PP)/(4.0d0*DX(2:SIZE_ST-1,PP))


VECT_R(1,1)=D_1(2,PP)+DF_1(2,PP)*U1(1,PP)*DT(PP)/(4.0d0*DX(2,PP))
VECT_R(1,2)=D_2(2,PP)+DF_2(2,PP)*U2(1,PP)*DT(PP)/(4.0d0*DX(2,PP))
VECT_R(1,3)=D_3(2,PP)+DF_3(2,PP)*U3(1,PP)*DT(PP)/(4.0d0*DX(2,PP))
VECT_R(1,4)=D_4(2,PP)+DF_4(2,PP)*U4(1,PP)*DT(PP)/(4.0d0*DX(2,PP))

VECT_R(SIZE_ST-2,1)=D_1(SIZE_ST-1,PP)-DF_1(SIZE_ST-1,PP)*U1(SIZE_ST,PP)*DT(PP)/(4.0d0*DX(SIZE_ST,PP))
VECT_R(SIZE_ST-2,2)=D_2(SIZE_ST-1,PP)-DF_2(SIZE_ST-1,PP)*U2(SIZE_ST,PP)*DT(PP)/(4.0d0*DX(SIZE_ST,PP))
VECT_R(SIZE_ST-2,3)=D_3(SIZE_ST-1,PP)-DF_3(SIZE_ST-1,PP)*U3(SIZE_ST,PP)*DT(PP)/(4.0d0*DX(SIZE_ST,PP))
VECT_R(SIZE_ST-2,4)=D_4(SIZE_ST-1,PP)-DF_4(SIZE_ST-1,PP)*U4(SIZE_ST,PP)*DT(PP)/(4.0d0*DX(SIZE_ST,PP))

VECT_R(2:SIZE_ST-3,1)=D_1(3:SIZE_ST-2,PP)
VECT_R(2:SIZE_ST-3,2)=D_2(3:SIZE_ST-2,PP)
VECT_R(2:SIZE_ST-3,3)=D_3(3:SIZE_ST-2,PP)
VECT_R(2:SIZE_ST-3,4)=D_4(3:SIZE_ST-2,PP)

DO CN=1,4

CALL TRIDIAG(DIAG_A(:,CN),DIAG_B(:,CN),DIAG_C(:,CN),VEC T_R(:,CN),SOL(:,CN),SIZE_ST-2,CODE)



END DO
DO j=2,SIZE_ST-1
!U1(J,PP)=SOL(SIZE_ST-2-(J-2),1)+VIS1(J,PP)
!U2(J,PP)=SOL(SIZE_ST-2-(J-2),2)+VIS2(J,PP)
!U3(J,PP)=SOL(SIZE_ST-2-(J-2),3)+VIS3(J,PP)
!U4(J,PP)=SOL(SIZE_ST-2-(J-2),4)+VIS4(J,PP)

U1(J,PP)=SOL(J-1,1)+VIS1(J,PP)
U2(J,PP)=SOL(J-1,2)+VIS2(J,PP)
U3(J,PP)=SOL(J-1,3)+VIS3(J,PP)
U4(J,PP)=SOL(J-1,4)+VIS4(J,PP)

END DO
CALL BOUNDARY_CONDITIONS

DO J=1,SIZE_ST

F1(j,PP)=U2(j,PP)

F2(j,PP)=U2(j,PP)**2/U1(j,PP)+GAMONE*(U4(j,PP)-(U1(j,PP)/2.0D0)*((U2(j,PP)/U1(j,PP))**2+(U3(j,PP)/U1(j,PP))**2))

F3(j,PP)=U2(j,PP)*U3(j,PP)/U1(j,PP)

F4(j,PP)=U2(j,PP)*U4(j,PP)/U1(j,PP)+U2(j,PP)*GAMONE*(U4(j,PP)/U1(j,PP)-0.50D0*((U2(j,PP)/U1(j,PP))**2+(U3(j,PP)/U1(j,PP))**2))

J2_STAR(j,PP)= GAMONE*(U4(j,PP)-(U1(j,PP)/2.0D0)*((U2(j,PP)/U1(j,PP))**2+(U3(j,PP)/U1(j,PP))**2))

J2_THETA(j,PP)=((U3(j,PP))**2)/U1(j,PP)*(1/MEAN_RADIUS(j,PP))

F_K(j,PP)=-U2(j,PP)*U3(j,PP)/U1(j,PP)*(1/MEAN_RADIUS(j,PP))

END DO

END SUBROUTINE CrankNicolson