Page 323 - Introduction to Computational Fluid Dynamics

P1: ICD
15:43
May 11, 2005
CB908/Date
0521853265appb
302
LW=0 0 521 85326 5 APPENDIX B. 1D CONDUCTION CODE
LE=0
IF(I.EQ.2)LW=1
IF(I.EQ.N-1)LE=1
DXE=X(I+1)-X(I)
DXEP=X(I+1)-XCF(I+1)
DXEM=XCF(I+1)-X(I)
DXW=X(I)-X(I-1)
DXWP=X(I)-XCF(I)
DXWM=XCF(I)-X(I-1)
C*** CALCULATE CELL FACE CONDUCTIVITY BY HARMONIC MEAN.
CONDSME=DXE/(DXEM/COND(I)+DXEP/COND(I+1))*(1-LE)+LE*COND(I+1)
CONDSMW=DXW/(DXWP/COND(I)+DXWM/COND(I-1))*(1-LW)+LW*COND(I-1)
AW(I)=CONDSMW*ACF(I)/DXW
AE(I)=CONDSME*ACF(I+1)/DXE
1 CONTINUE
RETURN
END
C *************************************************
SUBROUTINE BOUND
INCLUDE ’COM1D.FOR’
C *************************************************
STAB(2)=AW(2)
STAB(N-1)=AE(N-1)
C*** FOR I=1 BOUNDARY
IF(T1SPEC) THEN
SU(2)=SU(2)+AW(2)*(PSI*T(1)+(1-PSI)*(TO(1)-TO(2)))
SP(2)=SP(2)+AW(2)*PSI
AW(2)=0.0
ELSE IF(Q1SPEC) THEN
SU(2)=SU(2)+ACF(2)*(PSI*QB1+(1-PSI)*QB1O)
T(1)=QB1*ACF(2)/(AW(2)+SMALL)+T(2)
AW(2)=0.0
ELSE IF (H1SPEC) THEN
TERM1=HB1*ACF(2)+SMALL
TERM2=AW(2)+SMALL
TERM=1/(1/TERM1+ 1/TERM2)
SU(2)=SU(2)+PSI*TERM*TINF1+TERM1*(1-PSI)*(TINF1O-TO(1))
SP(2)=SP(2)+PSI*TERM
T(1)=(T(2)+TERM1/TERM2*TINF1)/(1+TERM1/TERM2)
AW(2)=0.0
ENDIF