Page 370 - Introduction to Computational Fluid Dynamics

P1: ICD/GKJ
0 521 85326 5
CB908/Date
0521853265appc
349
APPENDIX C. 2D CARTESIAN CODE
40 P2: IWV DO 41 J=2,JNM May 25, 2005 11:59
DO 41 I=2,INM
SU(I,J)=SU(I,J)+0.0
41 CONTINUE
GO TO 1000
C *** FOR KINETIC ENERGY
50 DO 51 J=2,JNM
DO 51 I=2,INM
IF(NTAG(I,J).EQ.1)GO TO 51
LW=NTAGW(I,J)
LE=NTAGE(I,J)
LS=NTAGS(I,J)
LN=NTAGN(I,J)
C EXCLUDE NEAR-WALL NODE
IF(LW.EQ.14.OR.LW.EQ.16)GO TO 51
IF(LE.EQ.24.OR.LE.EQ.26)GO TO 51
IF(LS.EQ.34.OR.LS.EQ.36)GO TO 51
IF(LN.EQ.44.OR.LN.EQ.46)GO TO 51
C PRODUCTION TERMS
DUDX=(FINTE(U,I,J)-FINTW(U,I,J))/DXP(I)
DVDX=(FINTE(V,I,J)-FINTW(V,I,J))/DXP(I)
DUDY=(FINTN(U,I,J)-FINTS(U,I,J))/DYP(J)
DVDY=(FINTN(V,I,J)-FINTS(V,I,J))/DYP(J)
TERM=2.0*(DUDX**2+DVDY**2)+(DUDY+DVDX)**2
IF(AXISYMM)TERM=TERM+2*(V(I,J)/R(J))**2
PROD(I,J)=VIST(I,J)*TERM
ENP=AMAX1(E(I,J),0.0)
SU(I,J)=PROD(I,J)*VOL(I,J) +SU(I,J)
SP(I,J)=RHO(I,J)**2*CMU*ENP/(VIST(I,J)+SMALL)*VOL(I,J)+SP(I,J)
51 CONTINUE
GO TO 1000
C *** FOR DISSIPATION
60 DO 61 J=2,JNM
DO 61 I=2,INM
IF(NTAG(I,J).EQ.1)GO TO 61
RHOP=RHO(I,J)
VOLP=VOL(I,J)
DPEP=ABS(D(I,J)/(E(I,J)+SMALL))
SU(I,J)=CD1*RHOP*DPEP*PROD(I,J)*VOLP +SU(I,J)
SP(I,J)=CD2*RHOP*DPEP*VOLP+SP(I,J)
61 CONTINUE
GO TO 1000