Page 324 - Introduction to Computational Fluid Dynamics

P1: ICD
0521853265appb
303
APPENDIX B. 1D CONDUCTION CODE
FOR I=N BOUNDARY
C*** CB908/Date 0 521 85326 5 May 11, 2005 15:43
IF(TNSPEC) THEN
SU(N-1)=SU(N-1)+AE(N-1)*(PSI*T(N)+(1-PSI)*(TO(N)-TO(N-1)))
SP(N-1)=SP(N-1)+AE(N-1)*PSI
AE(N-1)=0.0
ELSE IF(QNSPEC)THEN
SU(N-1)=SU(N-1)+ACF(N)*(PSI*QBN+(1-PSI)*QBNO)
T(N)=QBN*ACF(N)/(AE(N-1)+SMALL)+T(N-1)
AE(N-1)=0.0
ELSE IF(HNSPEC) THEN
TERM1=HBN*ACF(N)+SMALL
TERM2=AE(N-1)+SMALL
TERM=1/(1/TERM1+ 1/TERM2)
SU(N-1)=SU(N-1)+PSI*TERM*TINFN+TERM1*(1-PSI)*(TINFNO-TO(N))
SP(N-1)=SP(N-1)+PSI*TERM
T(N)=(T(N-1)+TERM1/TERM2*TINFN)/(1+TERM1/TERM2)
AE(N-1)=0.0
ENDIF
RETURN
END
C *************************************************
SUBROUTINE SOLVE
INCLUDE ’COM1D.FOR’
C *************************************************
DIMENSION AA(IT),BB(IT)
C*** ASSEMBLE SU AND SP TERMS
DO 1 I=2,N-1
IF(UNSTEADY)THEN
BP=RHO(I)*SPH(I)/DELT*VOL(I)
SP(I)=SP(I)+BP
SU(I)=SU(I)+(1-PSI)*(AE(I)*TO(I+1)+AW(I)*TO(I-1))
SU(I)=SU(I)+(BP-(1-PSI)*(AE(I)+AW(I)))*TO(I)
C CHECK FOR STABILITY CONDITION
TERM=BP-(1-PSI)*(AE(I)+AW(I)+STAB(I))
IF(TERM.LT.0.0)WRITE(*,*)’ COEF OF TPOLD IS NEGATIVE ATI=’,I
ENDIF
AP(I)=PSI*(AE(I)+AW(I))+SP(I)
C UNDER-RELAX
E=(1.-RP)/RP*AP(I)
AP(I)=AP(I)+E
SU(I)=SU(I)+E*T(I)
1 CONTINUE