CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Crank-Nicolson help request

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   November 3, 2012, 18:40
Default Crank-Nicolson help request
  #1
New Member
 
Join Date: Nov 2012
Posts: 2
Rep Power: 0
sezo85 is on a distinguished road
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
sezo85 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Saving ParaFoam views and case sail OpenFOAM Paraview & paraFoam 9 November 25, 2011 16:46
Crank Nicolson scheme msrinath80 OpenFOAM Running, Solving & CFD 6 November 14, 2010 14:59
Xwindows crash with paraview save srinath OpenFOAM Paraview & paraFoam 1 October 15, 2008 09:37
Crank Nicolson. ! Main CFD Forum 0 September 5, 2005 12:52
Convergence problem using crank nicolson method Carlos Gonzalez CFX 1 September 24, 2002 16:52


All times are GMT -4. The time now is 23:29.