|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？注册
x
**************************************************************************
* COND1DC (PART I) *
* *
* A FINITE VOLUME PROGRAM FOR THE SOLUTION OF ONE-DIMENSIONAL *
* STEADY CONDUCTION-TYPE PROBLEMS IN THE CARTESIAN COORDINATE SYSTEM *
* *
* DEVELOPED BY : B. RABI BALIGA *
* PREPARED FOR COMPUTATIONAL PHYSICS CLASS BY A. GHODS *
* *
* DATE : APRIL 21, 2000 *
* *
*------------------------------------------------------------------------*
* COND1DC IS A TWO-PART PROGRAM *
* *
* PART I: PROBLEM-INDEPENDENT PORTION OF COND1DC - IT CONSISTS *
* OF THE MAIN PROGRAM 'COND1DC' AND THE FOLLOWING *
* SUBROUTINES: *
* SETUP; COEFF; SOLVE; DETOUT; AND COEFFP *
**************************************************************************
* NOTATION :
* AN(J) - THE NORTH NEIGHBOR COEFFICIENT FOR THE J NODE
* AP(J) - COEFFICIENT FOR THE J NODE
* AS(J) - THE SOUTH NEIGHBOR COEFFICIENT FOR THE J NODE
* CON(J) - THE CONSTANT ASSOCIATED WITH THE DISCRETIZED
* EQUATION FOR THE J NODE
* CONVGE - LOGICAL EXPRESSION: DEFAULT VALUE IS .FALSE.:
* TERMINATES PROGRAM WHEN SET TO .TRUE.
* F(J,NF) - ANY ONE OF A MAXIMUM OF NV (NF = 1 TO NV) DEPENDENT
* VARIABLES OF INTEREST
* GAM(J) - DIFFUSION COEFFICIENT
* ITER - ITERATION COUNTER
* ITSRT - STARTING VALUE OF ITER IN A PARTICULAR RUN
* JREF - J INDEX OF A USER-SPECIFIED REFERENCE NODE
* LPRINT(NF) - LOGICAL VARIABLE: DEFAULT VALUE IS .FALSE.: IF SET
* TO .TRUE. A DETAILED PRINTOUT OF THE ASSOCIATED
* F(J,NF) IS OBTAINED
* LSOLVE(NF) - LOGICAL VARIABLE DEFAULT VALUE IS ,FALSE.: IF SET
* TO .TRUE., F(J,NF) IS SOLVED
* MAXIT - USER-SPECIFIED MAXIMUM NUMBER OF ITERATIONS DESIRED
* M1 - LAST NODE IN THE Y-DIR.
* M2 - SECOND TO LAST NODE IN THE Y-DIR.
* NF - REPRESENTS THE DEPENDENT VARIABLE BEING SOLVED
* NY - TOTAL # OF NODES IN THE Y DIRECTION: USER INPUT PARAMETER
* NV - TOTAL # OF DEPENDENT VARIABLES; USER INPUT PARAMETER
* PCOEFF - LOGICAL VARIABLE: DEFAULT VALUE IS .FALSE.: IF SET TO
* .TRUE., ALLOWS A DETAILED PRINTOUT OF COEFFICIENTS:
* USEFUL FOR DEBUGGING DURING PROGRAM DEVELOPMENT
* PT(J) - RECURRENCE VARIABLE USED IN TDMA ALGORITHM
* QT(J) - RECURRENCE VARIABLE USED IN TDMA ALGORITHM
* RELAX(NF) - DENOTES THE UNDER-RELAXATION PARAMETER FOR F(I,J,NF)
* SC(J) - CONSTANT THE LINEARIZED SOURCE TERM
* SP(J) - COEFFICIENT OF F(J,NF) IN THE LINEARIZED SOURCE TERM
* Y(J) - Y-COORDINATE OF JTH NODAL POINT
* YCV(J) - DISTANCE BETWEEN ADJACENT MAIN-GRID CV FACES
* IN THE Y-DIR.
* YDIF(J) - DISTANCE BETWEEN ADJACENT MAIN-GRID NODES (J - 1) & J
* YL - TOTAL LENGTH OF CALC. DOMAIN IN THE Y-DIR.
* YPOW - USED TO PRODUCE AN UNEVEN NODAL DISTRIBUTION IN Y-DIR.
* YV(J) - JTH MAIN-GRID CV FACE LOCATION
*****************************************************************************
* ---------------
PROGRAM COND1DC
* ---------------
* INCLUDING COMMON BLOCK
* ***********************
INCLUDE 'cblk3.for'
* DOMAIN DISCRETIZATION AND RELATED CALCULATIONS
* ----------------------------------------------
CALL GRID
CALL SETUP
* ASSIGN DEFAULT VALUES TO CONTROL VARIABLES & OTHER VARIABLES IF DESIRED
* -----------------------------------------------------------------------
CALL DEFVAL
* SOLUTION SEQUENCE FOR STEADY ONE-DIMENSIONAL CONDUCTION-TYPE PROBLEMS
* ---------------------------------------------------------------------
CALL START
10 ITER = ITER + 1
DO 20 NF = 1, NV
IF (LSOLVE(NF)) THEN
CALL GAMSOR
CALL COEFFS
CALL BOUND
CALL SOLVE
ENDIF
20 CONTINUE
CALL OUTPUT
IF ((.NOT.CONVGE).AND.(ITER.LT.MAXIT)) GO TO 10
STOP
END
**************************************************************************
* -----------------
SUBROUTINE SETUP
* -----------------
* THE FUNCTION OF THIS SUBROUTINE IS TO SET UP ALL GRID-RELATED QUANTITIES
* USED IN THE FORMULATION OF THE FINITE VOLUME METHOD
* -----------------------------------------------------------------------
* INCLUDING THE COMMON BLOCK
* **************************
INCLUDE 'cblk3.for'
* NOTATION :
* Y(J) - Y-COORDINATE OF JTH NODAL POINT
* YCV(J) - DISTANCE BETWEEN ADJACENT MAIN-GRID CV FACES
* IN THE Y-DIR.
* YDIF(J) - DISTANCE BETWEEN ADJACENT MAIN-GRID NODES (J - 1) & J
* YV(J) - JTH MAIN-GRID CV FACE LOCATION (USER SPECIFIED)
* --------------------------------------------------------------------
* THESE FOLLOWING VARIABLES ARE INITIALIZED, BUT NOT USED IN THE FVM
* ------------------------------------------------------------------
YDIF(1) = 0.0D+00
* SET UP GRID-RELATED QUANTITIES USING USER-SPECIFIED Y(J)'s
* NOTE: THIS FVM IS SETUP FOR TYPE-A GRID (REF. PATANKAR)
* -----------------------------------------------------------
YV(1) = Y(1)
DO 20 J = 2, M1
YV(J) = (Y(J-1) + Y(J)) * 0.5D+00
YDIF(J) = Y(J) - Y(J-1)
YCV(J-1) = YV(J) - YV(J-1)
20 CONTINUE
YV(M1+1) = Y(M1)
YCV(M1) = YV(M1+1) - YV(M1)
RETURN
END
********************************************************************************
* -----------------
SUBROUTINE COEFFS
* -----------------
* THE FUNCTION OF THIS SUBROUTINE IS TO EVALUATE THE COEFFICIENTS
* IN THE DISCRETIZED EQUATIONS
* ---------------------------------------------------------------
* INCLUDING COMMON BLOCK
* ***********************
INCLUDE 'cblk3.for'
* NOTATION :
* AN(J) - THE NORTH NEIGHBOR COEFFICIENT FOR THE J NODE
* AP(J) - COEFFICIENT FOR THE J NODE
* AS(J) - THE SOUTH NEIGHBOR COEFFICIENT FOR THE J NODE
* CON(J) - THE CONSTANT ASSOCIATED WITH THE DISCRETIZED
* EQUATION FOR THE J NODE
* --------------------------------------------------------------------
* THE SOUTH COEFFICIENT FOR NODE 1 IS ZERO
AS(1) = 0.0D+00
* THE NORTH COEFFICIENT FOR NODE M1 IS ZERO
AN(M1) = 0.0D+00
* CALCULATE THE OTHER NORTH AND SOUTH COEFFICIENTS
DO 10 J = 1, M2
JP1 = J + 1
* RESISTANCE-TYPE INTERPOLATION FOR THE DIFFUSION COEFFICIENTS
* NOTE: WITH THE TYPE-A GRID, THIS BECOMES THE HARMONIC MEAN
* - GAMN REFERS TO THE INTERPOLATED VALUE OF GAMMA AT THE
* NORTH FACE OF THE CV FOR NODE J
GAMN = 2.0D+00 * GAM(J) * GAM(JP1) / (GAM(J) + GAM(JP1)
+ + 1.0D-50)
AN(J) = GAMN / YDIF(JP1)
* NOTE: THE NORTH COEFFICIENT FOR NODE J EQUALS THE SOUTH COEFFICIENT
* FOR NODE J+1 - REMEMBER THAT THIS IS A CODUCTION CODE
AS(JP1) = AN(J)
10 CONTINUE
* CALCULATE THE PARTICULAR-POINT COEFFICIENTS AP(J) AND UNDER-RELAX:
* RELAX(NF) DENOTES THE UNDER-RELAXATION PARAMETER. ALSO CALCULATE
* THE CON(J) TERMS, AND THE UNDER-RELAXATION CONTRIBUTION
* NOTE: THE IMPLICIT UNDER-RELAXATION TECHNIQUE OF PATANKAR IS USED
REL = 1.0D+00 - RELAX(NF)
DO 20 J = 1, M1
AP(J) = (AN(J) + AS(J) - SP(J)*YCV(J)) / RELAX(NF)
CON(J) = SC(J) * YCV(J) + REL * AP(J) * F(J,NF)
20 CONTINUE
* IF AN OUTPUT OF THE COEFFS. IS DESIRED IN FILE COEFF.DAT,
* SET PCOEFF = .TRUE. IN SUBROUTINE START
IF (PCOEFF) CALL COEFFP
RETURN
END
********************************************************************************
* ----------------
SUBROUTINE SOLVE
* ----------------
* THIS SUBROUTINE IMPLEMENTS A TRI-DIAGONAL-MATRIX-ALGORITHM FOR THE
* SOLUTION OF THE LINEAR, OR NOMINALLY-LINEAR, DISCRETIZATION EQUATIONS
*--------------------------------------------------------------------------
* INCLUDING COMMON BLOCK
* ***********************
INCLUDE 'cblk3.for'
* PT AND QT ARE THE RECURRENCE VARIABLES FOR A TDMA
DIMENSION PT(NY), QT(NY)
* --------------------------
* STARTING VALUES FOR THE RECURRENCE VARIABLES
PT(1) = AN(1) / AP(1)
QT(1) = CON(1) / AP(1)
* CALCULATE OTHER RECURRENCE VARIABLES
DO 10 J = 2, M1
JM1 = J - 1
DENOM = AP(J) - AS(J) * PT(JM1)
PT(J) = AN(J) / DENOM
QT(J) = (CON(J) + AS(J) * QT(JM1)) / DENOM
10 CONTINUE
* INITIATE BACK-SUBSTITUTION.
F(M1,NF) = QT(M1)
* PERFORM THE BACK-SUBSTITUTION FOR J = M2, ............, 3, 1
DO 20 J = 1, M2
JJ = M1 - J
JJP1 = JJ + 1
F(JJ, NF) = PT(JJ) * F(JJP1, NF) + QT(JJ)
20 CONTINUE
RETURN
END
********************************************************************************
* -----------------
SUBROUTINE DETOUT
* -----------------
* THE FUNCTION OF THIS SUBROUTINE IS TO PROVIDE A DETAILED FIELD OUTPUT
* OF THE DEPENDENT VARIABLES OF INTEREST ( SET LPRINT(NF) = .TRUE. )
* ---------------------------------------------------------------------
* INCLUDING COMMON BLOCK
* ***********************
INCLUDE 'cblk3.for'
* FORMAT STATEMENTS
* -----------------
100 FORMAT( A4 , 4X, 6(I4, 8X) )
200 FORMAT( A7 , 6(1PD10.3, 2X) )
300 FORMAT( / A20, / A20 / )
400 FORMAT( /// 6(1H*), 3X, A5, I2, 2X, A6, I2, A2, 3X, 6(1H*) /
+ 9X, 19(1H-)//)
500 FORMAT( A4, 5X, 6(I4, 8X) )
600 FORMAT( 6X, 6(1PD10.3, 2X) )
700 FORMAT( 1X /)
* ------------------
WRITE(6,700)
WRITE(6,*) ' COMPUTATION OF STEADY CONDUCTION-TYPE PROBLEMS'
WRITE(6,*) ' IN A 1-D CARTESIAN COORDINATE SYSTEM'
WRITE(6,*) ' ----------------------------------------------'
JEND = 0
20 JBEG = JEND + 1
JEND = JEND + 6
IF (JEND .GT. M1) JEND = M1
WRITE(6,100) ' J =', (J, J = JBEG, JEND)
WRITE(6,200) ' Y(J) =', (Y(J), J = JBEG, JEND)
WRITE(6,200) ' YV(J)=', (YV(J), J = JBEG, JEND)
IF (JEND .LT. M1) GO TO 20
DO 50 NF = 1, NV
IF (LPRINT(NF)) THEN
WRITE(6,400) 'NF = ', NF, ' F(J, ', NF, ') '
JEND = 0
30 JBEG = JEND + 1
JEND = JEND + 6
IF (JEND.GT.M1) JEND = M1
WRITE(6,500) ' J =', (J, J= JBEG, JEND)
WRITE(6,600) (F(J,NF), J = JBEG, JEND)
IF (JEND.LT.M1) GO TO 30
ENDIF
50 CONTINUE
RETURN
END
********************************************************************************
* -----------------
SUBROUTINE COEFFP
* -----------------
* THE FUNCTION OF THIS SUBROUTINE IS TO OUTPUT THE COEFFICIENTS
* FOR THE DEPENDENT VARIABLES OF INTEREST, e.g. IF DEBUGGING IS REQUIRED
* NOTE: COEFFS. ARE OUTPUTTED TO FILE COEFF.DAT
* ----------------------------------------------------------------------
* INCLUDING COMMON BLOCK
* ***********************
INCLUDE 'cblk3.for'
* FORMAT STATEMENTS
* -----------------
100 FORMAT( A5, 1X, I2 // )
200 FORMAT( '1' /)
300 FORMAT( 2X, I2, 2X, 4(1PE10.3, 2X) )
400 FORMAT( A4, 2X, 4(2X, A4, 6X), / )
* ------------------------------------
OPEN( UNIT = 7, FILE = 'COEFF.DAT', STATUS = 'UNKNOWN')
REWIND 7
WRITE(7,200)
WRITE(7,100) '0NF =', NF
WRITE(7,400) ' J = ', ' AN ', ' AS ', ' AP ', ' CON '
DO 10 J = 1, M1
WRITE(7,300) J, AN(J), AS(J), AP(J), CON(J)
10 CONTINUE
CLOSE( UNIT = 7 )
RETURN
END
*******************************************************************************
**************************************************************************
* CD1DC (PART II) *
* *
* A FINITE VOLUME PROGRAM FOR THE SOLUTION OF ONE-DIMENSIONAL *
* STEADY CONDUCTION-TYPE PROBLEMS IN THE CARTESIAN COORDINATE SYSTEM *
* *
* DEVELOPED BY : B. RABI BALIGA *
* *
* PREPARED FOR STUDENTS IN THE GRADUATE COURSE ON COMPUTATIONAL PHYSICS *
* BY A. GHODS *
* DATE 21 APRIL 2000 *
* *
*------------------------------------------------------------------------*
* CD1DC IS A TWO-PART PROGRAM *
* *
* PART II: PROBLEM-DEPENDENT PORTION OF CD1DC - IT CONSISTS *
* OF THE FOLLOWING SUBROUTINES *
* GRID; DEFVAL; START; GAMSOR; BOUND; AND OUTPUT *
* *
* SPECIALIZATION: THIS VERSION IS DESIGNED TO SOLVE UNSTEADY *
* HEAT 1-D HEAT CONDUCTION IN THIN PLATE THERMAL *
* ISOLATED IN ALL OF ITS FACES EXCEP THOSE IN Y *
* DIRECTION *
**************************************************************************
* NOTATION :
* AN(J) - THE NORTH NEIGHBOR COEFFICIENT FOR THE J NODE
* AP(J) - COEFFICIENT FOR THE J NODE
* AS(J) - THE SOUTH NEIGHBOR COEFFICIENT FOR THE J NODE
* CON(J) - THE CONSTANT ASSOCIATED WITH THE DISCRETIZED
* EQUATION FOR THE J NODE
* CONVGE - LOGICAL EXPRESSION: DEFAULT VALUE IS .FALSE.:
* TERMINATES PROGRAM WHEN SET TO .TRUE.
* F(J,NF) - ANY ONE OF A MAXIMUM OF TEN (NF = 1 TO 10) DEPENDENT
* VARIABLES OF INTEREST
* GAM(J) - DIFFUSION COEFFICIENT
* ITER - ITERATION COUNTER
* ITSRT - STARTING VALUE OF ITER IN A PARTICULAR RUN
* JREF - J INDEX OF A USER-SPECIFIED REFERENCE NODE
* LPRINT(NF) - LOGICAL VARIABLE: DEFAULT VALUE IS .FALSE.: IF SET
* TO .TRUE. A DETAILED PRINTOUT OF THE ASSOCIATED
* F(J,NF) IS OBTAINED
* LSOLVE(NF) - LOGICAL VARIABLE DEFAULT VALUE IS ,FALSE.: IF SET
* TO .TRUE., F(J,NF) IS SOLVED
* MAXIT - USER-SPECIFIED MAXIMUM NUMBER OF ITERATIONS DESIRED
* M1 - LAST NODE IN THE Y-DIR.
* M2 - SECOND TO LAST NODE IN THE Y-DIR.
* NF - REPRESENTS THE DEPENDENT VARIABLE BEING SOLVED
* NY - TOTAL # OF NODES IN THE Y DIRECTION
* NV - TOTAL # OF DEPENDENT VARIABLES
* PCOEFF - LOGICAL VARIABLE: DEFAULT VALUE IS .FALSE.: IF SET TO
* .TRUE., ALLOWS A DETAILED PRINTOUT OF COEFFICIENTS:
* USEFUL FOR DEBUGGING DURING PROGRAM DEVELOPMENT
* PT(N) - RECURRENCE VARIABLE USED IN TDMA ALGORITHM
* QT(N) - RECURRENCE VARIABLE USED IN TDMA ALGORITHM
* RELAX(NF) - DENOTES THE UNDER-RELAXATION PARAMETER FOR F(I,J,NF)
* SC(J) - CONSTANT THE LINEARIZED SOURCE TERM
* SP(J) - COEFFICIENT OF F(J,NF) IN THE LINEARIZED SOURCE TERM
* Y(J) - Y-COORDINATE OF JTH NODAL POINT
* YCV(J) - DISTANCE BETWEEN ADJACENT MAIN-GRID CV FACES
* IN THE Y-DIR.
* YDIF(J) - DISTANCE BETWEEN ADJACENT MAIN-GRID NODES (J - 1) & J
* YL - TOTAL LENGTH OF CALC. DOMAIN IN THE Y-DIR.
* YPOW - USED TO PRODUCE AN UNEVEN NODAL DISTRIBUTION IN Y-DIR.
* YV(J) - JTH MAIN-GRID CV FACE LOCATION
*****************************************************************************
* ---------------
SUBROUTINE GRID
* ---------------
* THE PURPOSE OF THIS SUBROUTINE IS TO ENABLE THE USER TO INPUT THE
* COORDINATES OF THE GRID POINTS
* -----------------------------------------------------------------
* INCLUDING COMMON BLOCK
* ***********************
INCLUDE 'cblk3.for'
* NOTATION :
* M1 - LAST NODE IN THE Y-DIR.
* M2 - SECOND TO LAST NODE IN THE Y-DIR.
* YL - EXTENT FO THE CALCULATION DOMAIN IN THE Y-DIR.
* YPOW - EXPONENT IN A POWER-LAW TYPE DISTRIBUTION OF NODES
*
* NOTE: THIS IS ALSO A GOOD POINT TO ASSIGN A VALUE TO THE
* REFERENCE NODE, JREF
WRITE(*,*) ' M1 = ?'
WRITE(*,*) ' '
READ(*,*) M1
M2 = M1 - 1
JREF = M1
* INPUT THE LENGHT OF THE CONDUCTING MATERIAL, H, IN METERS
H = 100.D+00
* INPUT THE Y-EXTENT OF THE CALCULATION DOMAIN: HERE, WE'LL
YL = H
* GENERATE THE GRID: POWER-LAW TYPE GRID; YPOW = 1 IMPLIES UNIFORM
* GRID; YPOW > 1 CROWDS NODES IN THE VICINITY OF
* THE BOTTOM WALL
YPOW = 1.0D+00
Y(1) = 0.0D+00
Y(M1) = YL
DO 20 J = 2, M2
Y(J) = (DFLOAT(J-1) / DFLOAT(M1 - 1)) ** YPOW * YL
20 CONTINUE
RETURN
END
********************************************************************************
* -----------------
SUBROUTINE DEFVAL
* -----------------
* THE PURPOSE OF THIS SUBROUTINE IS TO ASSIGN DEFAULT VALUES TO CONTROL
* PARAMETERS AND ANY DESIRED VARIABLES
* ---------------------------------------------------------------------
* INCLUDING COMMON BLOCK
* ***********************
INCLUDE 'cblk3.for'
* NOTATION :
* CONVGE - TERMINATES EXECUTION OF THE PROGRAM WHEN SET TO TRUE
* ITER - ITERATION COUNTER
* ITSRT - STARTING VALUE OF ITER IN A PARTICULAR RUN
* LPRINT(NF) - LOGICAL VARIABLE: DEFAULT VALUE IS .FALSE.: IF SET
* TO .TRUE. A DETAILED PRINTOUT OF THE ASSOCIATED
* F(J,NF) IS OBTAINED
* LSOLVE(NF) - LOGICAL VARIABLE DEFAULT VALUE IS ,FALSE.: IF SET
* TO .TRUE., F(J,NF) IS SOLVED
* MAXIT - USER-SPECIFIED MAXIMUM NUMBER OF ITERATIONS DESIRED
* PCOEFF - LOGICAL VARIABLE: DEFAULT VALUE IS .FALSE.: IF SET TO
* .TRUE., ALLOWS A DETAILED PRINTOUT OF COEFFICIENTS:
* USEFUL FOR DEBUGGING DURING PROGRAM DEVELOPMENT
* RELAX(NF) - DENOTES THE UNDER-RELAXATION PARAMETER FOR F(I,J,NF)
* ---------------------------------------------------------------------
ITSRT= 0
ITER = ITSRT
MAXIT = 10
CONVGE = .FALSE.
PCOEFF = .TRUE.
DO 30 NF = 1, NV
LSOLVE(NF) = .FALSE.
LPRINT(NF) = .FALSE.
RELAX(NF) = 1.0D+00
* SET ALL F(J,NF) TO ZERO
DO 40 J = 1, M1
F(J,NF) = 0.0D+00
40 CONTINUE
30 CONTINUE
RETURN
END
*******************************************************************************
* ----------------
SUBROUTINE START
* ----------------
* THE PURPOSE OF THIS SUBROUTINE IS TO INPUT ALL PROBLEM-DEPENDENT
* PARAMETERS AND ASSIGN PROPER INITIAL OR GUESS VALUES TO F(J,NF)
* ----------------------------------------------------------------
* INCLUDING COMMON BLOCK
* ***********************
INCLUDE 'cblk3.for'
* ------------------------------------------------------------------
* SET LOGICAL PARAMETERS TO .TRUE. FOR NF = 1 : I.E. SOLVE FOR T(J)
LSOLVE(1) = .TRUE.
* SET THE RELAXATION PARAMETER FOR NF = 1
RELAX(1) = 1.0D+00
* ---------------------------------------
* READ IN VALUES OF VISC, DEN, AND DPDX
* HERE, VALUES CLOSE TO THOSE FOR AIR
* ARE SPECIFIED - THIS IS NOT GOING TO
* INFLUENCE THE NONDIMENSIONAL RESULTS
* - E.G. U/UAV AND FRE
* ---------------------------------------
* SPECIFY THERMAL CONDUCTIVITY IN ?
CONDUCT = 5.0
* SPECIFY THE VALUE OF HEAT GENERATION
HGEN=0
* SET GUESS VALUES FOR T(J) IN THE DOMAIN: THE WALL VALUE = 0.0
DO 10 J = 1, M1
T(J) = 0.0D+00
10 CONTINUE
* SPECIFY TEMPERATURE AT TOP AND BOTTOM BOUNDARY
TB=0.0D+00 ! B.C. AT THE BOTTOM
Tt=100.0D+00 ! B.C. AT THE TOP
T(1) =TB
T(M1) =TT
* NOTE: IF ISTART = 1, THEN THE PROBLEM IS RUN FROM SCRATCH (ITER = 0)
* IF ISTART > 1, THEN THE PROBLEM IS STARTED FROM THE RESULTS OF
* THE LAST RUN
WRITE(*,*) ' ISTART = ?'
WRITE(*,*) ' '
READ(*,*) ISTART
IF (ISTART.GT.1) THEN
OPEN (UNIT = 51, FILE = 'DATAIN.DAT', STATUS = 'UNKNOWN')
REWIND 51
READ(51,*) ITSRT
DO 20 J = 1, M1
READ(51,*) (F(J,NF), NF = 1, NV)
20 CONTINUE
ITER = ITSRT
WRITE(*,*) 'ITSRT = ', ITSRT, ' F(JREF,1) =', F(JREF,1)
CLOSE (UNIT = 51)
ENDIF
* SET MAXIT
WRITE(*,*) ' STARTING ITER = ', ITER
WRITE(*,*) ' MAXIT = ?'
WRITE(*,*) ' '
READ(*,*) MAXIT
RETURN
END
********************************************************************************
* -----------------
SUBROUTINE GAMSOR
* -----------------
* THE PURPOSE OF THIS SUBROUTINE IS TO ASSIGN ALL DIFFUSION COEFFICIENTS
* AND SOURCE TERMS: GAM(J), SC(J), AND SP(J) FOR EACH DEPENDENT VARIABLE
* THAT IS BEING SOLVED
* -----------------------------------------------------------------------
* INCLUDING COMMON BLOCK
* ***********************
INCLUDE 'cblk3.for'
* NOTE: HERE, ONLY NF = 1, I.E. T(J), IS BEING SOLVED
DO 30 J = 1, M1
GAM(J) = CONDUCT
SC(J) = HGEN
SP(J) = 0.0D+00
30 CONTINUE
RETURN
END
********************************************************************************
* ----------------
SUBROUTINE BOUND
* ----------------
* THE PURPOSE OF THIS SUBROUTINE IS TO INPUT BOUNDARY INFO., IF NEEDED
* --------------------------------------------------------------------
* INCLUDING COMMON BLOCK
* ***********************
INCLUDE 'cblk3.for'
* SET THE VALUE OF T AT THE BOTTOM WALL TO ZERO: T(1) = 0.0D+00
AP(1) = 1.0D+00
AN(1) = 0.0D+00
AS(1) = 0.0D+00
CON(1) = TB
* SET THE VALUE OF T AT THE TOP WALL TO Tt: T(J) = Tt
AP(M1) = 1.0D+00
AN(M1) = 0.0D+00
AS(M1) = 0.0D+00
CON(M1) =TT
* CALCULATE THE RESIDUE TO MONITOR THE CONVERGENCE IN OUTPUT
RESID = 0.0D+00
RES = 0.0D+00
DO 10 J = 1, M1
RES = AP(J) * T(J) - CON(J)
IF (J .NE. 1) THEN
RES = RES - AS(J) * T(J-1)
ENDIF
IF (J .NE. M1) THEN
RES = RES - AN(J) * T(J+1)
ENDIF
RESID = RESID + DABS(RES)
10 CONTINUE
RETURN
END
********************************************************************************
* -----------------
SUBROUTINE OUTPUT
* -----------------
* THIS SUBROUTINE IS MAINLY USED TO SET UP OUTPUT OF DESIRED RESULTS;
* BUT IT COULD ALSO BE USED TO SET CONVERGENCE CRITERIA; MONITOR
* CONVERGENCE; AND CALCULATE AUXILIARY RESULTS
* --------------------------------------------------------------------
* INCLUDING COMMON BLOCK
* ***********************
INCLUDE 'cblk3.for'
* OPEN THE MONIT.DAT FILE TO SAVE CONVERGENCE MONITORING INFO, IF NEEDED
* NOTE: THIS IS DONE ONLY FOR THE FIRST ITERATION OF EACH RUN
IF (ITER.LE.ITSRT+1) THEN
OPEN (UNIT = 11, FILE = 'MONIT.DAT', STATUS = 'UNKNOWN')
REWIND 11
ENDIF
WRITE(11,*) ' ITER =', ITER,' F(JREF,1) =', F(JREF,1)
WRITE(11,*)
* WRITE MONITORING INFO TO SCREEN IF DESIRED
WRITE(*,*) ' ITER =', ITER,' F(JREF,1) =', F(JREF,1)
WRITE(*,*)
* CONVERGENCE CRITERION
* NORMALIZE THE RESIDUE W.R.T. THE THERMAL. FLUX. AT THE WALL
WMFLUX = HGEN * YCV(1) + CONDUCT * (T(2) - T(1)) / YDIF(2)
RESID = RESID / WMFLUX
IF (RESID .LE. 1.0D-06) CONVGE = .TRUE.
IF ((ITER.EQ.MAXIT).OR.CONVGE) THEN
* CLOSE UNIT = 11 : THE MONIT FILE
CLOSE (UNIT = 11)
* OUTPUT ALL DESIRED DATA
* FIRST DUMP RESTART DATA ON TO DATAOUT.DAT
OPEN (UNIT = 50, FILE = 'DATAOUT.DAT', STATUS = 'UNKNOWN')
REWIND 50
WRITE(50,*) ITER
DO 20 J = 1,M1
WRITE(50,*) (F(J,INF), INF = 1, NV)
20 CONTINUE
CLOSE (UNIT = 50)
* PREPARE UNIT 6 FOR PLOTTING OUTPUT, IF DESIRED
OPEN (UNIT = 6, FILE = 'PLOTFILE.DAT', STATUS = 'UNKNOWN')
REWIND 6
WRITE(*,*)'Y(J) T(J)'
DO 500 J = 1, M1
WRITE(6,*) Y(J), T(J)
500 CONTINUE
CLOSE (UNIT = 6)
ENDIF
IF (CONVGE) THEN
WRITE(*,*) ' THE NUMERICAL SOLUTION CONVERGED.'
WRITE(*,*) ' NUMBER OF ITERATIONS : ', ITER
ELSE
CONTINUE
ENDIF
RETURN
END
*******************************************************************************
|
|
IP卡
狗仔卡
发表于 2003-4-12 03:44:53
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡
