- #1
Aerospark
- 2
- 0
Ok so basically I am trying to write code that can call from a fluid properties subroutine similar to refprop but it is called FLUID. I need to call the subroutine with temperature and pressure, and have it come back with density, compressibility, and etc. It should be able to do this, however I keep coming back with errors and I'm not sure why..
The error that I am getting according to the code description is due to being "out of range in G table."
Main Program:
The subroutine DRYAIR is one of many (theres one for O2,N2,CO2 etc.) But basically it fills G,TT,DD,GG,SAT,C, and M matrices up with coefficients which are used as transport equations in the fluid subroutine:
DRYAIR CODE:
Finally the code for the FLUID subroutine is here:
I'm really not sure why it is not returning me with density when I specify everything correctly, except I am not sure about how the common and implicit commands are affecting the code so it must be something to do with those.. The problem is that the G tables are "out of range"
Any help is greatly appreciated, and if i get it working I will be more than happy to send you a copy of it.
The error that I am getting according to the code description is due to being "out of range in G table."
Main Program:
Code:
PROGRAM GasProperties
INTEGER ENTRYdata,ERROR,IGOTO
LOGICAL VAPOR
DIMENSION PROPS(8)
COMMON/FLUIDC/ GAMMA,WL,WG,DENSL,DENSG,ENTL,ENTG,ENTHL,ENTHG
COMMON /FLUDPC/ G,TT,DD,GG,SAT,C,M,NG,N1,N2,N3,NS
REAL P,T,D,Z,S,H,CV,CP,GAM,Rgas,SONIC,MU,NU,K,Prandtl,RHO
CALL DRYAIR
WRITE(*,*) 'Enter P (Mpa):'
READ(*,*) P
WRITE(*,*) 'Enter T (Kelvin):'
READ(*,*) T
ENTRYdata=3
IGOTO=1
CALL FLUID(T,P,RHO,PROPS,8,ENTRYData,VAPOR,ERROR,IGOTO)
! PERFORM UNITS CONVERSIONS FOR OUTPUT AS FOLLOWS:
! Pressure from MPa to Pa
! Density from gm/cm^3 to kg/m^3
! Entropy from J/gmK to J/kgK
! Enthalpy from J/gmK to J/kgK
! Specific Heats from kJ/kgK to J/kgK
! Sonic velocity to m/s
! Viscosity (mu) from gm/cm-s to kg/m-s
! Thermal Conductivity from W/cmK to W/m-K
D = RHO*1000.0
Z=PROPS(1)
S = PROPS(2)*1000.0
H = PROPS(3)*1000.0
CV = PROPS(4)*1000.0
CP = PROPS(5)*1000.0
GAM = CP/CV
Rgas = CP-CV
SONIC = PROPS(6)*SQRT(1000.0)
MU = PROPS(7)*0.1
NU = MU/RHO
K = PROPS(8)*100.0
Prandtl = CP*MU/K
WRITE(*,*) 'Error Value: ', ERROR
WRITE(*,*)
WRITE(*,*) 'Density: ', RHO
WRITE(*,*) 'Compressibility Factor: ', Z
WRITE(*,*) 'Entropy ', S
WRITE(*,*) 'Enthalpy: ', H
WRITE(*,*) 'Specific Heat, Cv ', CV
WRITE(*,*) 'Specific Heat, Cp: ', CP
WRITE(*,*) 'Specific Heat Ratio, Gamma ', GAMMA
WRITE(*,*) 'Sonic Velocity: ', SONIC
WRITE(*,*) 'Dynamic Viscosity ', MU
WRITE(*,*) 'Kinematic Viscosity ', NU
WRITE(*,*) 'Thermal Conductivity ', K
WRITE(*,*) 'Prandtl Number, Pr ', Prandtl
END PROGRAM GasPropertiesThe subroutine DRYAIR is one of many (theres one for O2,N2,CO2 etc.) But basically it fills G,TT,DD,GG,SAT,C, and M matrices up with coefficients which are used as transport equations in the fluid subroutine:
DRYAIR CODE:
Code:
SUBROUTINE DRYAIR
INTERFACE
SUBROUTINE FLUID(T,P,RHO,PROPS,NP,ENTRYdata,VAPOR,ERROR,igoto)
REAL,INTENT(IN OUT):: t
REAL,INTENT(IN OUT):: p
REAL,INTENT(IN OUT):: rho
REAL,INTENT(OUT),DIMENSION(:):: props
INTEGER,INTENT(IN):: np
INTEGER,INTENT(IN):: ENTRYdata
LOGICAL,INTENT(OUT):: vapor
INTEGER,INTENT(OUT):: error
INTEGER,INTENT(IN):: igoto
END Subroutine Fluid
END INTERFACE
!**********************************************************************
REAL G(29,4),TT(27),DD(30),GG(27,30,8),SAT(40,8),C(8)
INTEGER M(4)
INTEGER ERROR,ENTRYdata
DIMENSION PROPS(8)
LOGICAL VAPOR
!**********************************************************************
DIMENSION X0001(100),X0101(100),X0201(100),X0301(100),X0401(100), &
& X0501(100),X0601(100),X0701(100),X0801(100),X0901(100), &
& X1001(100),X1101(100),X1201(100),X1301(100),X1401(100), &
& X1501(87)
DIMENSION Z(1587)
!**********************************************************************
COMMON /FLUDPC/ G,TT,DD,GG,SAT,C,M,NG,N1,N2,N3,NS
!**********************************************************************
EQUIVALENCE (X0001,Z(0001)), (X0101,Z(0101)), (X0201,Z(0201)), &
& (X0301,Z(0301)), (X0401,Z(0401)), (X0501,Z(0501)), &
& (X0601,Z(0601)), (X0701,Z(0701)), (X0801,Z(0801)), &
& (X0901,Z(0901)), (X1001,Z(1001)), (X1101,Z(1101)), &
& (X1201,Z(1201)), (X1301,Z(1301)), (X1401,Z(1401)), &
& (X1501,Z(1501))
!**********************************************************************
!**********************************************************************
DATA X0001/ 1.1300E+0, 1.2000E+0, 1.5000E+0, 1.9000E+0, 2.4000E+0,&
& 3.1000E+0, 3.9000E+0, 5.0000E+0, 6.3000E+0, 8.0000E+0, 1.0000E+1,&
& 1.1300E+1, 1.6892E+1, 1.7042E+1, 1.7598E+1, 1.8188E+1, 1.8778E+1,&
& 1.9422E+1, 2.0011E+1, 2.0670E+1, 2.1312E+1, 2.2013E+1, 2.2701E+1,&
& 2.3091E+1, 1.1897E+1, 1.2059E+1, 1.2751E+1, 1.3673E+1, 1.4830E+1,&
& 1.6462E+..... this goes on for a while, I left it out
!----------------------------------------------------------------------------
NG = 12
N1 = 12
N2 = 15
N3 = 8
NS = 8
DO 20 I=1,4
DO 10 J=1,12
G(J,I) = Z(J+(I-1)*12)
10 END DO
20 END DO
DO 30 I=1,12
TT(I) = Z(I+48)
30 END DO
DO 40 I=1,15
DD(I) = Z(I+60)
40 CONTINUE
KK=75
DO 70 I=1,8
DO 60 J=1,15
DO 50 K=1,12
KK=KK+1
GG(K,J,I) = Z(KK)
50 CONTINUE
60 CONTINUE
70 CONTINUE
DO 90 I=1,8
DO 80 J=1,8
SAT(J,I) = Z(J+(I-1)*8+1515)
80 CONTINUE
90 CONTINUE
DO 100 I=1,8
C(I) = Z(I+1579)
100 CONTINUE
IGOTO = 4
CALL FLUID(TEMP,PRES,DENS,PROPS,NP,ENTRYdata,VAPOR,ERROR,IGOTO)
RETURN
END SUBROUTINE DRYAIRFinally the code for the FLUID subroutine is here:
Code:
C WRITTEN BY: THEODORE E. FESSLER, NASA TM X-3572, AUGUST 1977
C REWRITTEN BY: LT. MARK D. KLEM & MARGARET P. PROCTOR, JANUARY 1986
C
C THE ORIGINAL PROGRAM WAS WRITTEN IN SFTRAN, BUT THIS PROGRAM IS
C REWRITTEN IN STANDARD FORTRAN 77 TO RUN PRIMARILY ON PERSONAL
C COMPUTERS AND CAN ALSO RUN ON MAINFRAMES. MUCH OF THE OLD CODE WAS
C COMMENTED OUT BUT LEFT INTACT FOR DOCUMENTATION AND REFERENCE
C PURPOSES.
C
SUBROUTINE FLUID (TEMP,PRES,DENS,PROPS,NP,ENTRY,VAPOR,ERROR,IGOTO)
C SUBROUTINE FLUID0 (G,NG,GG,TT,DD,N1,N2,N3,SAT,NS,C,L)
C...INTITIALIZING ENTRY FOR FLUID-PROPERTIES PACKAGE. SETS UP TABLE
C POINTERS FOR A PARTICULAR FLUID. SUBSEQUENT CALLS (TO FLUID) USE
C THESE POINTERS.
C...CALLING ARGUMENTS:
C G = (NON-DIMENSIONAL) TABLES IN TEMPERATURE ONLY.
C G(1,1) = ARGUMENT VALUES (TEMPERATURE).
C G(1,2) = ENTROPY, S0.
C G(1,3) = ENERGY, U0.
C G(1,4) = SPECIFIC HEAT AT CONSTANT VOLUME, CV0.
C NG = LENGTH OF EACH G TABLE.
C GG = TABLES IN TEMPERATURE AND DENSITY.
C GG(1,1,1) = Z1, THE MODEL-DEPARTURE FACTOR.
C GG(1,1,2) = Z2 FACTOR FOR ENTROPY (S).
C GG(1,1,3) = Z3 FACTOR FOR ENTHALPY (H).
C GG(1,1,4) = Z4 FACTOR FOR SPECIFIC HEAT (CV).
C GG(1,1,5) = Z5 FACTOR FOR SPECIFIC HEAT (CP).
C GG(1,1,6) = Z6 FACTOR FOR VELOCITY OF SOUND.
C GG(1,1,N) = OTHER FACTORS (EG. TRANSPORT PROPERTIES).
C TT = TEMPERATURE ARGUMENTS (REDUCED) FOR GG TABLES.
C DD = DENSITY ARGUMENTS (REDUCED) FOR GG TABLES.
C N1 = LENGTH OF TT TABLE AND FIRST DIMENSION OF GG TABLE.
C N2 = LENGTH OF DD TABLE AND SECOND DIMENSION OF GG TABLE.
C N3 = THIRD DIMENSION (NUMBER OF) GG TABLES.
C SAT = (NON-DIMENSIONAL) SATURATION TABLES.
C SAT(1,1) = SATURATION TEMPERATURE.
C SAT(1,2) = SATURATION PRESSURE.
C SAT(1,3) = SATURATED LIQUID DENSITY.
C SAT(1,4) = SATURATED GAS DENSITY.
C SAT(1,5) = SATURATED LIQUID ENTROPY.
C SAT(1,6) = SATURATED GAS ENTROPY.
C SAT(1,7) = SATURATED LIQUID ENTHALPY.
C SAT(1,8) = SATURATED GAS ENTHALPY.
C NS = LENGTH OF EACH SATURATION TABLE.
C C = MISCELLANEOUS CONSTANTS:
C C(1) = SPECIFIC GAS CONSTANT, R.
C C(2) = CRITICAL TEMPERATURE, TC.
C C(3) = CRITICAL PRESSURE, PC.
C C(4) = CRITICAL DENSITY, DC (OR PSEUDO VALUE FOR BEST RANGE).
C C(5) = CONVERGENCE REQUIREMENT FOR TEMPERATURE.
C C(6) = CONVERGENCE REQUIREMENT FOR DENSITY.
C C(7) = CONVERGENCE REQUIREMENT FOR ENTROPY.
C C(8) = CONVERGENCE REQUIREMENT FOR ENTHALPY.
C L = MEMORY (4 WORDS) FOR TABLE INDICES.
C...VARIABLES IN COMMON BLOCK /FLUIDC/:
C GAMMA = RATIO OF SPECIFIC HEATS, CP/CV.
C WL = MASS-FRACTION IN THE LIQUID PHASE.
C WG = MASS-FRACTION IN THE GAS PHASE.
C DENSL = DENSITY OF SATURATED LIQUID, SAME UNITS AS DC.
C DENSG = DENSITY OF SATURATED GAS, SAME UNITS AS DC.
C ENTL = ENTROPY OF SATURATED LIQUID, SAME UNITS AS R.
C ENTG = ENTROPY OF SATURATED GAS, SAME UNITS AS R.
C ENTHL = ENTHALPY OF SATURATED LIQUID, SAME UNITS AS R*TC.
C ENTHG = ENTHALPY OF SATURATED GAS, SAME UNITS AS R*TC.
C...EXTERNAL SUBROUTINES:
C NTRP; CALCULATION OF INTERPOLATION COEFFICIENTS FOR SINGLE-
C VARIABLE INTERPOLATIONS,
C NTRP1; CALCULATION OF INTERPOLATION COEFFICIENTS FOR FIRST
C VARIABLE IN 2-VARIABLE INTERPOLATIONS,
C NTRP2; CALCULATION OF INTERPOLATION COEFFICIENTS FOR SECOND
C VARIABLE IN 2-VARIABLE INTERPOLATIONS,
C GNTRP; GETS VALUE IN SINGLE-VARIABLE INTERPOLATIONS,
C GGNTRP; GETS VALUE IN 2-VARIABLE INTERPOLATIONS,
C CUBIC; SOLVES VAN DER WAALS' EQUATION OF STATE FOR DENSITY.
C INCLUDE (TYPE STATEMENTS, ETC. FOR INITIALIZING ENTRY)
C INCLUDE (TYPE STATEMENTS, ETC. FOR MAIN ENTRY)
C DEFINITION (TYPE STATEMENTS, ETC. FOR MAIN ENTRY)
C DEFINITION (TYPE STATEMENTS, ETC. FOR INITIALIZING ENTRY)
REAL PROPS(NP)
INTEGER ENTRY,ERROR,ERROUT,BIT1,BIT2,BIT3,BIT4,BIT5
INTEGER DUMN0,DUMIA
DIMENSION DUMF(35),DUMFF(27,30),DUMA(1)
LOGICAL VAPOR,GAS
INTEGER L(4)
C REAL G(NG,4),GG(N1,N2,N3),TT(N1),DD(N2),SAT(NS,8),C(8),ROOTS(2)
REAL G(29,4),GG(27,30,8),TT(27),DD(30),SAT(40,8),C(8),ROOTS(2)
COMMON /FLUDPC/ G,TT,DD,GG,SAT,C,L,NG,N1,N2,N3,NS
COMMON /FLUIDC/ GAMMA,WL,WG,DENSL,DENSG,ENTL,ENTG,ENTHL,ENTHG
DATA BIT1,BIT2,BIT3,BIT4,BIT5,MAX2,MAX3/1,2,4,8,16,12,12/,
1 MAX4,MAX5,DXMAX,HUGE,SMALL,BIG/12,12,3.0,1.0E37,.001,1000./,
2 ERROUT / 6 /
C
C INCLUDE (FUNCTION DEFINITIONS)
RLIMIT(XMIN,X,XMAX)=AMAX1(XMIN,AMIN1(X,XMAX))
TF(P,D)=(P+3.0*D**2)*(3.0-D)/(8.0*Z*D)
PF(T,D)=8.0*Z*D*T/(3.0-D)-3.0*D**2
HF(P,D)=1.125*(P/(3.0*D)-D)
SF(D)=ALOG((3.0-D)/D)
CPF(T,D)=1.0/(1.0-D*(3.0-D)**2/(4.0*T))
A2F(T,D)=0.375*GAMMA*R*TC*(8.0*T/(3.0-D)*(1.0+D/(3.0-D))-6.0*D)
C...INITIALIZE NORMALIZING VALUES.
IF (IGOTO.EQ.4) THEN
R=C(1)
TC=C(2)
PC=C(3)
DC=C(4)
C RETURN
C
ELSE IF (IGOTO.EQ.1) THEN
C ENTRY FLUID (TEMP,PRES,DENS,PROPS,NP,ENTRY,VAPOR,ERROR)
C...MAIN ENTRY FOR PROPERTY CALCULATIONS.
C...CALLING ARGUMENTS:
C TEMP = FLUID TEMPERATURE, SAME UNITS AS TC.
C PRES = FLUID PRESSURE, SAME UNITS AS PC.
C DENS = FLUID DENSITY, SAME UNITS AS DC.
C PROPS = OTHER FLUID PROPERTIES:
C PROPS(1) = COMPRESSIBILITY, PRES/(R*DENS*TEMP).
C PROPS(2) = ENTROPY, SAME UNITS AS R.
C PROPS(3) = ENTHALPY, SAME UNITS AS R*TC.
C PROPS(4) = SPECIFIC HEAT, CV, SAME UNITS AS R.
C PROPS(5) = SPECIFIC HEAT, CP, SAME UNITS AS R.
C PROPS(6) = SONIC VELOCITY, SAME UNITS AS SQRT(R*TC).
C PROPS(7) = VISCOSITY, GRAMS/CM-SEC
C PROPS(8) = THERMAL CONDUCTIVITY, WATTS/CM-K
C NP = NUMBER OF PROPERTIES TO BE CALCULATED.
C ENTRY = INTEGER THAT SPECIFIES WHICH VARIABLES ARE INPUT:
C = 1 IF TEMPERATURE AND DENSITY ARE GIVEN.
C = 2 IF PRESSURE AND DENSITY ARE GIVEN.
C = 3 IF TEMPERATURE AND PRESSURE ARE GIVEN.
C = 4 IF PRESSURE AND ENTROPY ARE GIVEN.
C = 5 IF PRESSURE AND ENTHALPY ARE GIVEN.
C VAPOR = .TRUE. IF THE FLUID IS SATURATED. IN THAT CASE,
C VALUES OF LIQUID AND GAS PHASES ARE GIVEN IN THE
C COMMON BLOCK /FLUIDC/.
C ERROR = ERROR FLAGS (BITS -- LEAST SIGNIFICANT = 1):
C ** IF ERROR = 0, ALL IS WELL **
C BIT 1 = OUT OF RANGE IN SAT TABLE.
C BIT 2 = OUT OF RANGE IN G TABLE.
C BIT 3 = OUT OF RANGE IN TT TABLE.
C BIT 4 = OUT OF RANGE IN DD TABLE.
C BIT 5 = CONVERGENCE NOT ACHIEVED IN SOLVING INVERSE FUNCTION.
C PROCEDURE (MAIN ENTRY INITIALIZATION)
IF (NP.GT.N3) THEN
WRITE (ERROUT,1) NP,N3
STOP
ENDIF
IF (ENTRY.LT.1 .OR. ENTRY.GT.5) THEN
WRITE (ERROUT,2)
STOP
ENDIF
KS=0
K0=0
K1=0
K2=0
ERROR=0
ARG0=HUGE
ARG1=HUGE
ARG2=HUGE
T=TEMP/TC
D=DENS/DC
P=PRES/PC
C
C DO CASE (ENTRY,5)
IF (ENTRY.EQ.1) THEN
C CASE 1
C...GIVEN TEMPERATURE AND DENSITY.
IF (T.LT.1.0) THEN
C PROCEDURE (TEST FOR VAPOR CONDITION AT (T,D))
C PROCEDURE (SET UP FOR INTERPOLATION OF SATURATION DATA AT T)
CALL NTRP (SAT(1,1),T,NS,L(1),KS,1,DUMF,DUMFF,DUMVAL)
C PROCEDURE (GET SATURATED LIQUID (DL) AND GAS (DG) DENSITIES)
C CALL GNTRP (4,SAT(1,3),DL)
CALL NTRP (DUMA,DUMX,DUMN0,DUMIA,DUMOUT,4,SAT(1,3),
& DUMFF,DL)
C CALL GNTRP (4,SAT(1,4),DG)
CALL NTRP (DUMA,DUMX,DUMN0,DUMIA,DUMOUT,4,SAT(1,4),
& DUMFF,DG)
VAPOR=D.GT.DG .AND. D.LT.DL
ELSE
VAPOR=.FALSE.
ENDIF
IF (VAPOR) THEN
C PROCEDURE (GET SATURATION PRESSURE, PSAT)
C CALL GNTRP (4,SAT(1,2),PSAT)
CALL NTRP (DUMA,DUMX,DUMN0,DUMIA,DUMOUT,4,
& SAT(1,2),DUMFF,PSAT)
P=PSAT
ELSE
C PROCEDURE (GET Z BY DOUBLE INTERPOLATION AT (T,D))
C PROCEDURE (SET UP FOR DOUBLE INTERPOLATION AT (T,D))
IF (ARG1.NE.T) THEN
C CALL NTRP1 (TT,T,N1,L(3),K1,2)
CALL NTRP (TT,T,N1,L(3),K1,2,DUMF,DUMFF,DUMVAL)
ARG1=T
ENDIF
IF (ARG2.NE.D) THEN
C CALL NTRP2 (DD,D,N2,L(4),K2,3)
CALL NTRP (DD,D,N2,L(4),K2,3,DUMF,DUMFF,DUMVAL)
ARG2=D
ENDIF
C CALL GGNTRP (5,GG(1,1,1),Z)
CALL NTRP (DUMA,DUMX,DUMN0,DUMIA,DUMOUT,5,DUMF,
& GG(1,1,1),Z)
Z=RLIMIT(SMALL,Z,BIG)
P=PF(T,D)
ENDIF
PRES=P*PC
ELSE IF (ENTRY.EQ.2) THEN
C CASE 2
C...GIVEN PRESSURE AND DENSITY.
IF (P.LT.1.0) THEN
C PROCEDURE (TEST FOR VAPOR CONDITION AT (P,D))
C PROCEDURE (SET UP FOR INTERPOLATION OF SATURATION DATA AT P)
CALL NTRP (SAT(1,2),AMAX1(P,SAT(1,2)),NS,L(1),KS,1,
& DUMF,DUMFF,DUMVAL)
C PROCEDURE (GET SATURATED LIQUID (DL) AND GAS (DG) DENSITIES)
C CALL GNTRP (4,SAT(1,3),DL)
CALL NTRP (DUMA,DUMX,DUMN0,DUMIA,DUMOUT,4,
& SAT(1,3),DUMFF,DL)
C CALL GNTRP (4,SAT(1,4),DG)
CALL NTRP (DUMA,DUMX,DUMN0,DUMIA,DUMOUT,4,
& SAT(1,4),DUMFF,DG)
VAPOR=D.GT.DG .AND. D.LT.DL
ELSE
VAPOR=.FALSE.
ENDIF.... keeps going with more cases but that is the jist of itI'm really not sure why it is not returning me with density when I specify everything correctly, except I am not sure about how the common and implicit commands are affecting the code so it must be something to do with those.. The problem is that the G tables are "out of range"
Any help is greatly appreciated, and if i get it working I will be more than happy to send you a copy of it.
