💻 Fortran Source Code Library
Run calculations locally on your own machine. View code structure, read technical explanations, and download compilation packages including sample input files.
Steam Tables (NBS/NRC Helmholtz Solver)
Core Numerical Engine in Fortran 90 • 2 total downloads
steam.f90 / steam_calc.f90
! =========================================================================
! Source File: steam.f90
! =========================================================================
subroutine base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
!*****************************************************************************80
!
!! base() calculates quantities associated with the base Helmholtz function.
!
! Discussion:
!
! The equation for the base Helmholtz function AB(T,RHO) is:
!
! AB(T,RHO) = R * T * (
! - ln ( 1 - y )
! - ( beta - 1 ) / ( 1 - y )
! + ( alpha + beta + 1 ) / ( 2 * ( 1 - y )^2 )
! + 4 * y * ( ( Bbar / b ) - gamma )
! - 0.5 * ( alpha - beta + 3 )
! + ln ( RHO * R * T / P0 ) )
! (Equation 2)
! where
!
! y = b * rho / 4,
! alpha = 11,
! beta = 133/3,
! gamma = 7/2,
! P0 = 0.101325 MegaPascals = 1 atm
!
! and
!
! b(T) = b1 * ln(T/T0) + sum(j=0,1,3,5) b(j)*(T0/T)^j (Equation 3)
!
! Bbar(T) = sum(j=0,1,2,4) B(j)*(T0/T)^j (Equation 4).
!
! where
!
! T0=647.073 K and the coefficients b(j) and B(j) are
!
! j b(j) B(j)
! -- ----------- ----------
! 0 0.7478629 1.1278334
! 1 -0.3540782 -0.5944001
! 2 0 -5.010996
! 3 0.007159876 0
! 4 0 0.63684256
! 5 -0.003528426 0
!
! For the derived quantities, the following relations are used:
!
! Pressure: PB = RHO^2 * dAB/dRHO
! Density derivative: DPDRB = 2*PB/RHO + RHO^2 * d2AB/dRHO2
! Temperature derivative: DPDTB = RHO^2 * d2AB/(dRHO dT)
! Specific entropy: SB = ( UB - AB ) / T
! Specific internal energy: UB = AB + T * SB
! Specific enthalpy: HB = UB + PB / RHO
! Specific heat capacity
! at constant volume: CVB = - T * d2AB/dT2
! Specific Gibbs function: GB = AB + PB / RHO
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 03 February 2002
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Input:
!
! real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! real ( kind = rk8 ) RHO, the density, in G/CM3.
!
! Output:
!
! real ( kind = rk8 ) AB, the base value of the Helmholtz function,
! in KJ/kg.
!
! real ( kind = rk8 ) CVB, the base value of the isochoric (constant
! volume) heat capacity, in KJ/(kg degrees Kelvin).
!
! real ( kind = rk8 ) DPDRB, the base value of the partial
! derivative dP(T,RHO)/dRHO, with T held fixed, in (MegaPascals CM3)/G.
!
! real ( kind = rk8 ) DPDTB, the base value of the partial
! derivative dP(T,RHO)/dT, with RHO held fixed, in
! MegaPascals/degrees Kelvin.
!
! real ( kind = rk8 ) GB, the base value of the Gibbs free energy,
! in KJ/kg.
!
! real ( kind = rk8 ) HB, the base value of enthalpy, in KJ/kg.
!
! real ( kind = rk8 ) PB, the base pressure, in MegaPascals.
!
! real ( kind = rk8 ) SB, the base value of entropy,
! in KJ/(kg degrees Kelvin).
!
! real ( kind = rk8 ) UB, the base value of internal energy, in KJ/kg.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) ab
real ( kind = rk8 ), parameter :: alpha = 11.0D+00
real ( kind = rk8 ) b1
real ( kind = rk8 ) b1t
real ( kind = rk8 ) b1tt
real ( kind = rk8 ) b2
real ( kind = rk8 ) b2t
real ( kind = rk8 ) b2tt
real ( kind = rk8 ), parameter :: beta = 44.333333333333D+00
real ( kind = rk8 ) cvb
real ( kind = rk8 ) dpdrb
real ( kind = rk8 ) dpdtb
real ( kind = rk8 ) dz
real ( kind = rk8 ) dz0
real ( kind = rk8 ), parameter :: gamma = 3.5D+00
real ( kind = rk8 ) gascon
real ( kind = rk8 ) gb
real ( kind = rk8 ) hb
real ( kind = rk8 ), parameter :: p_zero = 0.101325D+00
real ( kind = rk8 ) pb
real ( kind = rk8 ) rho
real ( kind = rk8 ) sb
real ( kind = rk8 ) t
real ( kind = rk8 ) ub
real ( kind = rk8 ) x
real ( kind = rk8 ) y
real ( kind = rk8 ) z
real ( kind = rk8 ) z0
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BASE - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BASE - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was RHO = ', rho
stop
end if
!
! Compute auxilliary quantities for Equation 2.
!
call bb ( t, b1, b2, b1t, b2t, b1tt, b2tt )
y = 0.25D+00 * b1 * rho
x = 1.0D+00 - y
!
! Evaluate Equation 2.
!
ab = - log ( 1.0D+00 - y ) &
- ( beta - 1.0D+00 ) / ( 1.0D+00 - y ) &
+ ( alpha + beta + 1.0D+00 ) / ( 2.0D+00 * ( 1.0D+00 - y )**2 ) &
+ 4.0D+00 * y * ( ( b2 / b1 ) - gamma ) &
- 0.5D+00 * ( alpha - beta + 3.0D+00 ) &
+ log ( rho * gascon() * t / p_zero )
!
! Determine quantities defined in terms of AB.
!
pb = ( 1.0D+00 + alpha * y + beta * y**2 ) / ( 1.0D+00 - y )**3 &
+ 4.0D+00 * y * ( b2 / b1 - gamma )
z0 = ( 1.0D+00 + alpha * y + beta * y**2 ) / ( 1.0D+00 - y )**3
z = z0 + 4.0D+00 * y * ( b2 / b1 - gamma )
dz0 = ( alpha + 2.0D+00 * beta * y ) / ( 1.0D+00 - y )**3 &
+ 3.0D+00 * ( 1.0D+00 + alpha * y + beta * y**2 ) / ( 1.0D+00 - y )**4
dz = dz0 + 4.0D+00 * ( b2 / b1 - gamma )
gb = ab + pb
ub = - t * b1t * ( pb - 1.0D+00 - rho * b2 ) / b1 - rho * t * b2t
hb = pb + ub
!
! An incorrect version of this equation began:
!
! cvb = 2.0D+00 * ub + ( pb - 1.0D+00 ) &
!
! and caused me no end of trouble. My fault, JVB, 03 February 2002
!
cvb = 2.0D+00 * ub + ( z0 - 1.0D+00 ) &
* ( ( t * b1t / b1 )**2 - t**2 * b1tt / b1 ) &
- rho * t**2 * ( b2tt - gamma * b1tt ) - ( t * b1t / b1 )**2 * y * dz0
dpdtb = pb / t + rho * ( 0.25D+00 * ( dz0 + 4.0D+00 * ( b2 / b1 - gamma ) ) &
* b1t + b2t - b2 / b1 * b1t )
sb = ub - ab
dpdrb = pb + y * ( dz0 + 4.0D+00 * ( b2 / b1 - gamma ) )
!
! Assign dimensions.
!
ab = gascon() * t * ab
cvb = gascon() * cvb
dpdrb = gascon() * t * dpdrb
dpdtb = gascon() * t * rho * dpdtb
gb = gascon() * t * gb
hb = gascon() * t * hb
pb = gascon() * t * rho * pb
sb = gascon() * sb
ub = gascon() * t * ub
return
end
subroutine bb ( t, b1, b2, b1t, b2t, b1tt, b2tt )
!*****************************************************************************80
!
!! bb() calculates the B's of equations 3 and 4.
!
! Discussion:
!
! Here
!
! b(T) = b1 * ln(T/T0) + sum(j=0,1,3,5) b(j)*(T0/T)**j (Equation 3)
!
! Bbar(T) = sum(j=0,1,2,4) B(j)*(T0/T)**j (Equation 4).
!
! where
!
! T0 = 647.073 K
!
! and the coefficients b(j) and B(j) are
!
! j b(j) B(j)
! -- ----------- ----------
! 0 0.7478629 1.1278334
! 1 -0.3540782 -0.5944001
! 2 0 -5.010996
! 3 0.007159876 0
! 4 0 0.63684256
! 5 -0.003528426 0
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 13 January 2013
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) B1, the coefficient b from equation 3,
! in CM3/G.
!
! Output, real ( kind = rk8 ) B2, the coefficient Bbar from equation 4,
! in CM3/G.
!
! Output, real ( kind = rk8 ) B1T, the derivative dB1/dT,
! in (CM3)/(G Degrees Kelvin).
!
! Output, real ( kind = rk8 ) B2T, the derivative dB2/dT,
! in (CM3)/(G Degrees Kelvin).
!
! Output, real ( kind = rk8 ) B1TT, the second derivative of B1 with
! respect to T, in (CM3)/(G (Degrees Kelvin)^2 ).
!
! Output, real ( kind = rk8 ) B2TT, the second derivative of B2 with
! respect to T, in (CM3)/(G (Degrees Kelvin)^2 ).
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) b1
real ( kind = rk8 ) b1t
real ( kind = rk8 ) b1tt
real ( kind = rk8 ) b2
real ( kind = rk8 ) b2t
real ( kind = rk8 ) b2tt
real ( kind = rk8 ), parameter, dimension ( 10 ) :: bp = (/ &
0.7478629D+00, -0.3540782D+00, 0.0D+00, 0.0D+00, &
0.007159876D+00, 0.0D+00, -0.003528426D+00, 0.0D+00, &
0.0D+00, 0.0D+00 /)
real ( kind = rk8 ), parameter, dimension ( 10 ) :: bq = (/ &
1.1278334D+00, 0.0D+00, -0.5944001D+00, -5.010996D+00, &
0.0D+00, 0.63684256D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00 /)
integer i
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_ref = 647.073D+00
real ( kind = rk8 ) v(10)
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BB - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Set V(I) = ( T_REF / T )**(I-1).
!
v(1) = 1.0D+00
do i = 2, 10
v(i) = v(i-1) * t_ref / t
end do
!
! Set B1, B1T, B1TT.
!
b1 = bp(1) + bp(2) * log ( 1.0D+00 / v(2) )
b1t = bp(2) * v(2) / t_ref
b1tt = 0.0D+00
do i = 3, 10
b1 = b1 + bp(i) * v(i-1)
b1t = b1t - dble ( i - 2 ) * bp(i) * v(i-1) / t
b1tt = b1tt + bp(i) * dble ( i - 2 )**2 * v(i-1) / t**2
end do
b1tt = b1tt - ( b1t / t )
!
! Set B2, B2T, B2TT.
!
b2 = bq(1)
b2t = 0.0D+00
b2tt = 0.0D+00
do i = 3, 10
b2 = b2 + bq(i) * v(i-1)
b2t = b2t - dble ( i - 2 ) * bq(i) * v(i-1) / t
b2tt = b2tt + bq(i) * dble ( i - 2 )**2 * v(i-1) / t**2
end do
b2tt = b2tt - ( b2t / t )
return
end
subroutine corr ( t, p, p_consistent, rhol, rhov, delg )
!*****************************************************************************80
!
!! corr() evaluates an adjustment to the Gibbs function.
!
! Discussion:
!
! CORR is given T and P at or near the vapor pressure and evaluates
! the corresponding liquid and vapor densities, and the residual
! function DELG = (GL-GV)/(R*T) where GL and GV are the Gibbs functions
! for the liquid and vapor phases, respectively.
!
! These quantities are used to calculate a correction to the vapor
! pressure or the vapor temperature.
!
! The states corresponding to the coexisting phases of liquid
! and vapor for the temperature range from the triple point
! to within 0.5 C of the critical point 0.01 <= t <= tk-0.5 C
! have been determined in exact accord with the Gibbs condition
! of phase equilibrium: DELG = G(g)-G(l) = 0, P, t constant,
! where G(g) and G(l) are the values of the Gibbs function
! for saturated gas and liquid respectively.
!
! For the region (tk-t)<=0.5 C, an exact solution for the
! Helmholtz function yields values of density for the saturated
! liquid that are shifted to lower values. Also, the isotherms
! in the pressure-density plane and the Gibbs function-density
! plane are nearly flat, so that it is difficult to obtain
! solutions. As an alternative to exact solution, the power
! law equation is used to define states:
!
! rho(gas) = 0.322 - 0.657 * (1 - T/647.126)**0.325 (g/cm3).
! rho(liq) = 0.322 + 0.657 * (1 - T/647.126)**0.325 (g/cm3).
!
! In a poor instance of programming, the input pressure was
! originally overwritten on output by a value consistent with
! the computed densities. This causes no end of misunderstandings,
! since other routines expect the value of pressure to be input
! only. The code is now revised so that there is an input P
! and an output P. In a huff, JVB 05 February 2002.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 05 February 2002
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the vapor temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) P, the vapor pressure, in MegaPascals.
!
! Output, real ( kind = rk8 ) P_CONSISTENT, the vapor pressure, in MegaPascals,
! consistent with RHOL and RHOV. This is equal to the input value of
! P unless 646.3 <= T.
!
! Input/output, real ( kind = rk8 ) RHOL, the liquid density, in G/CM3.
! On input, if RHOL is positive, it is used as an initial
! estimate for the iteration.
!
! Input/output, real ( kind = rk8 ) RHOV, the vapor density, in G/CM3.
! On input, if RHOV is positive, it is used as an initial
! estimate for the iteration.
!
! Output, real ( kind = rk8 ) DELG, the residual function (GL-GV)/(R*T),
! where GL is the liquid Gibbs function, GV the vapor Gibbs function,
! dimensionless. If 646.3 < T, DELG is 0.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) a
real ( kind = rk8 ) ab
real ( kind = rk8 ) ar
real ( kind = rk8 ) cd
real ( kind = rk8 ) cjth
real ( kind = rk8 ) cjtt
real ( kind = rk8 ) cp
real ( kind = rk8 ) cv
real ( kind = rk8 ) cvb
real ( kind = rk8 ) cvr
logical, parameter :: debug = .false.
real ( kind = rk8 ) delg
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdrb
real ( kind = rk8 ) dpdrr
real ( kind = rk8 ) dpdt
real ( kind = rk8 ) dpdtb
real ( kind = rk8 ) dpdtr
real ( kind = rk8 ) g
real ( kind = rk8 ) gascon
real ( kind = rk8 ) gb
real ( kind = rk8 ) gl
real ( kind = rk8 ) gr
real ( kind = rk8 ) gv
real ( kind = rk8 ) h
real ( kind = rk8 ) hb
real ( kind = rk8 ) hr
real ( kind = rk8 ) p
real ( kind = rk8 ) p_consistent
real ( kind = rk8 ), parameter :: p_crit = 22.055D+00
real ( kind = rk8 ) pb
real ( kind = rk8 ) pr
real ( kind = rk8 ) rho
real ( kind = rk8 ) rhol
real ( kind = rk8 ), parameter :: rho_min = 1.0D-08
real ( kind = rk8 ) rhov
real ( kind = rk8 ) rho_start
real ( kind = rk8 ) s
real ( kind = rk8 ) sb
real ( kind = rk8 ) sr
real ( kind = rk8 ) t
real ( kind = rk8 ) tau
real ( kind = rk8 ), parameter :: t_crit = 647.1260000001D+00
real ( kind = rk8 ) u
real ( kind = rk8 ) ub
real ( kind = rk8 ) ur
p_consistent = p
!
! Initialize output quantities.
!
delg = 0.0D+00
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'CORR - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative pressures.
!
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'CORR - Fatal error!'
write ( *, '(a)' ) ' The input pressure P must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was P = ', p
stop
end if
if ( t <= 646.3D+00 ) then
if ( rhol <= 0.0D+00 ) then
rho_start = 1.11D+00 - 0.0004D+00 * t
else
rho_start = rhol
end if
call dense ( p_consistent, t, rho_start, rho, dpdr )
call therm ( t, rho, a, cjth, cjtt, cd, cv, dpdr, dpdt, g, h, &
p_consistent, s, u )
rhol = rho
gl = g
if ( rhov <= 0.0D+00 ) then
rho_start = p_consistent / ( gascon() * t )
else
rho_start = rhov
end if
call dense ( p_consistent, t, rho_start, rho, dpdr )
rho = max ( rho, rho_min )
call therm ( t, rho, a, cjth, cjtt, cp, cv, dpdr, dpdt, g, h, &
p_consistent, s, u )
rhov = rho
gv = g
delg = ( gl - gv ) / ( gascon() * t )
p_consistent = p
if ( debug ) then
write ( *, '(a,g14.6)' ) ' CORR - RHOL = ', rhol
write ( *, '(a,g14.6)' ) ' RHOV = ', rhov
end if
else if ( t <= t_crit ) then
if ( debug ) then
write ( *, '(a)' ) 'CORR - Twilight zone'
end if
delg = 0.0D+00
tau = 0.657128D+00 * ( 1.0D+00 - t / t_crit )**0.325D+00
rhol = 0.322D+00 + tau
rhov = 0.322D+00 - tau
rho = rhov
call base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
call resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
p_consistent = pb + pr
else
if ( debug ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'CORR - Weirdo zone'
end if
rhol = 0.322D+00
rhov = 0.322D+00
p_consistent = p_crit
delg = 0.0D+00
end if
return
end
subroutine cp_values ( n, tc, p, cp )
!*****************************************************************************80
!
!! cp_values() returns some values of the specific heat at constant pressure.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, pages 229-237.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) CP, the specific heat at constant pressure,
! in KJ/(kg K).
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 24
real ( kind = rk8 ) cp
real ( kind = rk8 ), save, dimension ( n_data ) :: cpvec = (/ &
4.228D+00, 2.042D+00, 1.975D+00, 2.013D+00, 2.040D+00, &
2.070D+00, 2.135D+00, 2.203D+00, 2.378D+00, 2.541D+00, &
2.792D+00, 2.931D+00, 4.226D+00, 4.223D+00, 4.202D+00, &
4.177D+00, 4.130D+00, 4.089D+00, 4.053D+00, 4.021D+00, &
3.909D+00, 3.844D+00, 3.786D+00, 2.89D+00 /)
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 5.0D+00, 10.0D+00, 50.0D+00, &
100.0D+00, 200.0D+00, 300.0D+00, 400.0D+00, 500.0D+00, &
1000.0D+00, 1500.0D+00, 2000.0D+00, 5000.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 100.0D+00, 200.0D+00, 300.0D+00, 350.0D+00, &
400.0D+00, 500.0D+00, 600.0D+00, 850.0D+00, 1100.0D+00, &
1600.0D+00, 2000.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
cp = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
cp = cpvec(n)
end if
return
end
subroutine dense ( p, t, rho_start, rho, dpdr )
!*****************************************************************************80
!
!! dense() computes the density for a given pressure and temperature.
!
! Discussion:
!
! The use of the variable RHO_START for two opposing purposes is
! poor practice and will be corrected one of these days. Meanwhile,
! the algorithm's behavior, particularly in the two-phase region,
! is very suspect.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 19 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) P, the pressure, in MegaPascals.
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO_START, an initial guess for the density,
! in G/CM3. The value of RHO_START also signals whether a vapor or liquid
! calculation is to be done. If DPDR is computed negative, then for
! 0.2967 <= RHO_START, liquid is assumed, otherwise gas.
!
! Output, real ( kind = rk8 ) RHO, the density for the given
! pressure and temperature, in G/CM3.
!
! Output, real ( kind = rk8 ) DPDR, the partial derivative
! dP(T,RHO)/dRHO, with T held fixed, in (MegaPascals CM3)/G.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) ab
real ( kind = rk8 ) ar
real ( kind = rk8 ) cvb
real ( kind = rk8 ) cvr
real ( kind = rk8 ) dp
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdrb
real ( kind = rk8 ) dpdrr
real ( kind = rk8 ) dpdtb
real ( kind = rk8 ) dpdtr
real ( kind = rk8 ) dpdx
real ( kind = rk8 ) errtol
real ( kind = rk8 ) gb
real ( kind = rk8 ) gr
real ( kind = rk8 ) hb
real ( kind = rk8 ) hr
integer it
integer, parameter :: it_max = 50
real ( kind = rk8 ) p
real ( kind = rk8 ) pb
real ( kind = rk8 ) pp
real ( kind = rk8 ) pr
real ( kind = rk8 ) rho
real ( kind = rk8 ), parameter :: rho_max = 1.9D+00
real ( kind = rk8 ), parameter :: rho_min = 1.0D-08
real ( kind = rk8 ) rho_start
real ( kind = rk8 ) sb
real ( kind = rk8 ) sr
real ( kind = rk8 ) t
real ( kind = rk8 ) ub
real ( kind = rk8 ) ur
real ( kind = rk8 ) x
errtol = sqrt ( epsilon ( errtol ) )
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'DENSE - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative pressures.
!
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'DENSE - Fatal error!'
write ( *, '(a)' ) ' The input pressure P must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was P = ', p
stop
end if
rho = rho_start
rho = max ( rho, rho_min )
rho = min ( rho, rho_max )
do it = 1, it_max
call resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
call base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
pp = pb + pr
dpdr = dpdrb + dpdrr
!
! Check for negative DP/DRho, which characterizes the two-phase region.
!
if ( dpdr <= 0.0D+00 ) then
if ( 0.2967D+00 <= rho_start ) then
rho = rho * 1.02D+00
else
rho = rho * 0.98D+00
end if
if ( it <= 10 ) then
cycle
end if
end if
dpdx = 1.1D+00 * dpdr
dpdx = max ( dpdx, 0.01D+00 )
dp = abs ( 1.0D+00 - pp / p )
if ( dp <= errtol .or. &
( 0.3D+00 < rho .and. dp <= errtol ) .or. &
( 0.7D+00 < rho .and. dp <= 10.0D+00 * errtol ) ) then
return
end if
x = ( p - pp ) / dpdx
if ( 0.1D+00 < abs ( x ) ) then
x = x * 0.1D+00 / abs ( x )
end if
rho = rho + x
rho = max ( rho, rho_min )
rho = min ( rho, rho_max )
end do
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'DENSE - Warning!'
write ( *, '(a)' ) ' The iteration did not converge.'
write ( *, '(a,i6)' ) ' Number of iterations was ', it_max
write ( *, '(a,g14.6)' ) ' Last iterate was ', rho
return
end
subroutine dielectric ( t, rho, eps )
!*****************************************************************************80
!
!! dielectric() returns the static dielectric constant.
!
! Discussion:
!
! According to the IAPS, the equation used is valid in the range
!
! 273.15 degrees Kelvin <= T <= 823.15 degrees K
! 0 MegaPascals <= P <= 500 MegaPascals.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, page 266.
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO, the density, in G/CM3.
!
! Output, real ( kind = rk8 ) EPS, the dielectric constant, dimensionless.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: npol_c = 4
real ( kind = rk8 ), parameter, dimension ( 10 ) :: a = (/ &
7.62571D+00, 244.003D+00, -140.569D+00, 27.7841D+00, -96.2805D+00, &
41.7909D+00, -10.2099D+00, -45.2059D+00, 84.6395D+00, -35.8644D+00 /)
real ( kind = rk8 ) c(0:npol_c)
real ( kind = rk8 ) eps
real ( kind = rk8 ) rho
real ( kind = rk8 ) t
real ( kind = rk8 ) t_copy
real ( kind = rk8 ), parameter :: t_max = 823.15D+00
real ( kind = rk8 ), parameter :: t_min = 273.15D+00
real ( kind = rk8 ), parameter :: t_ref = 298.15D+00
real ( kind = rk8 ) t_star
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'DIELECTRIC - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'DIELECTRIC - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was RHO = ', rho
stop
end if
t_copy = t
t_copy = min ( t_copy, t_max )
t_copy = max ( t_copy, t_min )
t_star = t_copy / t_ref
c(0) = 1.0D+00
c(1) = a(1) / t_star
c(2) = ( a(2) / t_star ) + a(3) + a(4) * t_star
c(3) = ( a(5) / t_star ) + a(6) * t_star + a(7) * t_star**2
c(4) = ( a(8) / t_star**2 ) + ( a(9) / t_star ) + a(10)
call r8poly_val_horner ( npol_c, c, rho, eps )
return
end
subroutine dielectric_values ( n, tc, p, eps )
!*****************************************************************************80
!
!! dielectric_values() returns some values of the static dielectric constant.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 03 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, page 266.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) EPS, the dielectric constant, dimensionless.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 15
real ( kind = rk8 ) eps
real ( kind = rk8 ), save, dimension ( n_data ) :: epsvec = (/ &
88.29D+00, 90.07D+00, 92.02D+00, 95.14D+00, 100.77D+00, &
78.85D+00, 70.27D+00, 62.60D+00, 55.78D+00, 44.31D+00, &
35.11D+00, 20.40D+00, 1.17D+00, 1.11D+00, 1.08D+00 /)
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
100.0D+00, 500.0D+00, 1000.0D+00, 2000.0D+00, 5000.0D+00, &
100.0D+00, 100.0D+00, 100.0D+00, 100.0D+00, 100.0D+00, &
100.0D+00, 100.0D+00, 100.0D+00, 100.0D+00, 100.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
25.0D+00, 50.0D+00, 75.0D+00, 100.0D+00, 150.0D+00, &
200.0D+00, 300.0D+00, 400.0D+00, 500.0D+00, 600.0D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
eps = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
eps = epsvec(n)
end if
return
end
function gascon ( )
!*****************************************************************************80
!
!! gascon() returns the value of the specific gas constant.
!
! Discussion:
!
! The specific gas constant R is related to the universal gas
! constant R-bar = 8.31441 J/(mol degrees Kelvin) by the molar mass
! M = 18.0152 g/mol:
!
! R = R-bar / M.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 13 January 2013
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Output, real ( kind = rk8 ) GASCON, the value of the specific gas
! constant, in J/(g degrees Kelvin).
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) gascon
gascon = 0.461522D+00
return
end
subroutine ideal ( t, ai, cpi, cvi, gi, hi, si, ui )
!*****************************************************************************80
!
!! ideal() computes ideal gas thermodynamic properties of water.
!
! Discussion:
!
! Values for thermodynamic properties of water in the ideal
! gas state were reported by Woolley. The formula for the ideal gas
! term of the Helmholtz function approximates a term by term summation of
! contributions from each of the rotation and vibration states.
! The formula, equation #6 in the reference, is:
!
! A(ideal)(T) = -R * T * ( 1 + ( C(1)/Tr + C(2) ) * ln(Tr)
! + Sum ( 3 <= I <= 18) C(I) * Tr**(I-6)
!
! where Tr=T/100 K. The C(i) are tabulated coefficients. Equation
! 6 can be used for temperatures below 3000 K, and is accurate to
! within the tolerance of the gas constant for 50<=T<=2000 K.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 13 January 2013
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) AI, the Helmholtz function, in KJ/kg.
!
! Output, real ( kind = rk8 ) CPI, the heat capacity at constant pressure,
! in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) CVI, the heat capacity at constant volume,
! in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) GI, the Gibbs free energy, in KJ/kg.
!
! Output, real ( kind = rk8 ) HI, the enthalpy, in KJ/kg.
!
! Output, real ( kind = rk8 ) SI, the entropy, in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) UI, the internal energy, in KJ/kg.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) ai
real ( kind = rk8 ), parameter, dimension ( 18 ) :: c = (/ &
19.730271018D+00, 20.9662681977D+00, -0.483429455355D+00, &
6.05743189245D+00, 22.56023885D+00, -9.87532442D+00, &
-4.3135538513D+00, 0.458155781D+00, -0.047754901883D+00, &
0.0041238460633D+00, -0.00027929052852D+00, 0.14481695261D-04, &
-0.56473658748D-06, 0.16200446D-07, -0.3303822796D-09, &
0.451916067368D-11, -0.370734122708D-13, 0.137546068238D-15 /)
real ( kind = rk8 ) cpi
real ( kind = rk8 ) cvi
real ( kind = rk8 ) gascon
real ( kind = rk8 ) gi
real ( kind = rk8 ) hi
integer i
real ( kind = rk8 ) si
real ( kind = rk8 ) t
real ( kind = rk8 ) tt
real ( kind = rk8 ) ui
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'IDEAL - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
tt = t / 100.0D+00
gi = - ( c(1) / tt + c(2) ) * log ( tt )
do i = 3, 18
gi = gi - c(i) * tt**(i-6)
end do
hi = c(2) + c(1) * ( 1.0D+00 - log ( tt ) ) / tt
do i = 3, 18
hi = hi + dble ( i - 6 ) * c(i) * tt**(i-6)
end do
cpi = c(2) - c(1) / tt
do i = 3, 18
cpi = cpi + dble ( ( i - 6 ) * ( i - 5 ) ) * c(i) * tt**(i-6)
end do
ai = gi - 1.0D+00
ui = hi - 1.0D+00
cvi = cpi - 1.0D+00
si = hi - gi
!
! Assign dimensions.
!
ai = gascon() * t * ai
cpi = gascon() * cpi
cvi = gascon() * cvi
gi = gascon() * t * gi
hi = gascon() * t * hi
si = gascon() * si
ui = gascon() * t * ui
return
end
subroutine prandtl ( t, p, pr )
!*****************************************************************************80
!
!! prandtl() computes the Prandtl number.
!
! Discussion:
!
! This routine was NOT working properly for large pressures,
! because the routine CORR was changing the input value of P.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 17 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) P, the pressure, in MegaPascals.
!
! Output, real ( kind = rk8 ) PR, the Prandtl number, dimensionless.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) a
real ( kind = rk8 ) cjth
real ( kind = rk8 ) cjtt
real ( kind = rk8 ) cp
real ( kind = rk8 ) cv
logical, parameter :: debug = .false.
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdt
real ( kind = rk8 ) eta
real ( kind = rk8 ) g
real ( kind = rk8 ) h
real ( kind = rk8 ) lambda
real ( kind = rk8 ) p
real ( kind = rk8 ) pr
real ( kind = rk8 ) rho
real ( kind = rk8 ) rhol
real ( kind = rk8 ) rho_start
real ( kind = rk8 ) rhov
real ( kind = rk8 ) s
real ( kind = rk8 ) t
real ( kind = rk8 ) t2
real ( kind = rk8 ) u
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PRANDTL - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative pressures.
!
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PRANDTL - Fatal error!'
write ( *, '(a)' ) ' The input pressure P must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was P = ', p
stop
end if
!
! Compute the density.
!
if ( debug ) then
write ( *, * ) 'PRANDTL - Call TSAT, with P = ', p
end if
rhol = 0.0D+00
rhov = 0.0D+00
call tsat ( p, t2, rhol, rhov )
if ( debug ) then
write ( *, * ) 'PRANDTL - T2 = ', t2
end if
if ( t < t2 ) then
rho_start = 1.9D+00
else
rho_start = 0.01D+00
end if
call dense ( p, t, rho_start, rho, dpdr )
if ( debug ) then
write ( *, * ) 'PRANDTL - RHO = ', rho
end if
!
! Now from T and RHO, compute CP, ETA and LAMBDA.
!
call therm ( t, rho, a, cjth, cjtt, cp, cv, dpdr, dpdt, g, h, p, s, u )
call viscosity ( t, rho, eta )
call thercon ( t, rho, lambda )
if ( debug ) then
write ( *, '(7f10.4)' ) t, p, rho, eta, cp, lambda, pr
end if
pr = eta * cp / lambda
return
end
subroutine prandtl_values ( n, tc, p, pr )
!*****************************************************************************80
!
!! prandtl_values() returns some values of the Prandtl number for testing.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, page 265.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) PR, the Prandtl number, dimensionless.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 35
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ) pr
real ( kind = rk8 ), save, dimension ( n_data ) :: prvec = (/ &
13.50D+00, 13.48D+00, 13.46D+00, 13.39D+00, 13.27D+00, &
13.15D+00, 13.04D+00, 12.93D+00, 12.83D+00, 12.73D+00, &
12.63D+00, 12.53D+00, 12.43D+00, 12.34D+00, 12.25D+00, &
12.08D+00, 11.92D+00, 11.77D+00, 11.62D+00, 11.48D+00, &
11.36D+00, 11.23D+00, 11.12D+00, 10.91D+00, 10.72D+00, &
10.55D+00, 6.137D+00, 3.555D+00, 2.378D+00, 1.000D+00, &
0.974D+00, 0.960D+00, 0.924D+00, 0.899D+00, 0.882D+00 /)
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
1.0D+00, 5.0D+00, 10.0D+00, 25.0D+00, 50.0D+00, &
75.0D+00, 100.0D+00, 125.0D+00, 150.0D+00, 175.0D+00, &
200.0D+00, 225.0D+00, 250.0D+00, 275.0D+00, 300.0D+00, &
350.0D+00, 400.0D+00, 450.0D+00, 500.0D+00, 550.0D+00, &
600.0D+00, 650.0D+00, 700.0D+00, 800.0D+00, 900.0D+00, &
1000.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 25.0D+00, 50.0D+00, 75.0D+00, 100.0D+00, &
150.0D+00, 200.0D+00, 400.0D+00, 600.0D+00, 800.0D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
pr = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
pr = prvec(n)
end if
return
end
subroutine psat ( t, p, rhol, rhov )
!*****************************************************************************80
!
!! psat() calculates the vapor pressure, and the liquid and vapor densities.
!
! Discussion:
!
! These quantities correspond to the input temperature T, corrected
! so that the Gibbs functions for liquid and vapor phase are
! equal to within a tolerance.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the vapor temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) P, the vapor pressure, in MegaPascals.
!
! Output, real ( kind = rk8 ) RHOL, the liquid density, in G/CM3.
!
! Output, real ( kind = rk8 ) RHOV, the vapor density, in G/CM3.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) bot
real ( kind = rk8 ) delg
real ( kind = rk8 ) dp
real ( kind = rk8 ) errtol
real ( kind = rk8 ) gascon
integer it
integer, parameter :: it_max = 100
real ( kind = rk8 ) p
real ( kind = rk8 ) p_consistent
real ( kind = rk8 ) p_old
real ( kind = rk8 ) rhol
real ( kind = rk8 ) rhov
real ( kind = rk8 ) t
!
! Ensure that output quantities are initialized,, obliterating any
! input values.
!
p = 0.0D+00
rhol = 0.0D+00
rhov = 0.0D+00
!
! Set the error tolerance.
!
errtol = 100.0D+00 * sqrt ( epsilon ( errtol ) )
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PSAT - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Get an estimate for the saturation pressure.
!
call psat_est ( t, p )
dp = 0.0D+00
do it = 1, it_max
call corr ( t, p, p_consistent, rhol, rhov, delg )
bot = ( rhol - rhov ) / ( rhol * rhov )
if ( abs ( bot ) < errtol ) then
write ( *, * ) 'PSAT - Warning, what is this?'
bot = sign ( errtol, bot )
end if
dp = delg * gascon() * t / bot
p_old = p
p = p + dp
if ( abs ( dp ) <= errtol * ( abs ( p ) + 1.0D+00 ) ) then
return
end if
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PSAT - Warning!'
write ( *, '(a)' ) ' The iterates have become nonpositive.'
write ( *, '(a,i6)' ) ' Iteration number = ', it
write ( *, '(a,g14.6)' ) ' Last iterate was ', p
write ( *, '(a,g14.6)' ) ' Previous iterate was ', p_old
write ( *, '(a,g14.6)' ) ' Last correction was ', dp
write ( *, '(a)' ) ' Trying to recover...'
p = 0.5D+00 * p_old
end if
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PSAT - Fatal error!'
write ( *, '(a)' ) ' The iterates have become nonpositive.'
write ( *, '(a,i6)' ) ' Iteration number = ', it
write ( *, '(a,g14.6)' ) ' Last iterate was ', p
write ( *, '(a,g14.6)' ) ' Previous iterate was ', p_old
write ( *, '(a,g14.6)' ) ' Last correction was ', dp
stop
end if
end do
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PSAT - Warning!'
write ( *, '(a)' ) ' The iteration did not converge.'
write ( *, '(a,i6)' ) ' The number of iterations was ', it_max
write ( *, '(a,g14.6)' ) ' Convergence tolerance was ', errtol
write ( *, '(a,g14.6)' ) ' Last iterate was ', p
write ( *, '(a,g14.6)' ) ' Last correction was ', dp
return
end
subroutine psat_est ( t, p )
!*****************************************************************************80
!
!! psat_est() makes a rough estimate of the vapor pressure.
!
! Discussion:
!
! The calculation agrees with tabulated data to within
! 0.02% for temperature to within a degree or so of the critical
! temperature. The approximate vapor pressure can be refined
! by imposing the condition that the Gibbs functions of the vapor
! and liquid phases be equal.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 21 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) P, the vapor pressure, in MegaPascals.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ), parameter, dimension ( 8 ) :: a = (/ &
-7.8889166D+00, 2.5514255D+00, -6.716169D+00, 33.239495D+00, &
-105.38479D+00, 174.35319D+00, -148.39348D+00, 48.631602D+00 /)
real ( kind = rk8 ) b
integer i
real ( kind = rk8 ) p
real ( kind = rk8 ) q
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_ref = 647.25D+00
real ( kind = rk8 ) v
real ( kind = rk8 ) w
real ( kind = rk8 ) z
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PSAT_EST - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
if ( t <= 314.0D+00 ) then
p = 0.1D+00 * exp ( 6.3573118D+00 - 8858.843D+00 / t &
+ 607.56335D+00 * t**( -0.6D+00 ) )
else
v = t / t_ref
w = abs ( 1.0D+00 - v )
b = 0.0D+00
do i = 1, 8
z = i
b = b + a(i) * w**( ( z + 1.0D+00 ) / 2.0D+00 )
end do
q = b / v
p = 22.093D+00 * exp ( q )
end if
return
end
subroutine psat_values ( n, tc, p )
!*****************************************************************************80
!
!! psat_values() returns some values of the saturation pressure.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, pages 9-15.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the saturation pressure, in bar.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 12
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
0.0061173D+00, 0.0065716D+00, 0.0087260D+00, 0.12344D+00, 1.0132D+00, &
2.3201D+00, 4.7572D+00, 15.537D+00, 39.737D+00, 85.838D+00, &
165.21D+00, 220.55D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.01D+00, 1.0D+00, 5.0D+00, 50.0D+00, 100.0D+00, &
125.0D+00, 150.0D+00, 200.0D+00, 250.0D+00, 300.0D+00, &
350.0D+00, 373.976D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
end if
return
end
subroutine r8poly_val_horner ( n, c, x, cx )
!*****************************************************************************80
!
!! r8poly_val_horner() evaluates a polynomial using Horner's method.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 08 August 1999
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the dimension of C.
!
! Input, real ( kind = rk8 ) C(0:N), the polynomial coefficients.
! C(I) is the coefficient of X^I.
!
! Input, real ( kind = rk8 ) X, the point at which the polynomial
! is to be evaluated.
!
! Output, real ( kind = rk8 ) CX, the value of the polynomial at X.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer n
real ( kind = rk8 ) c(0:n)
real ( kind = rk8 ) cx
integer i
real ( kind = rk8 ) x
cx = c(n)
do i = n - 1, 0, -1
cx = cx * x + c(i)
end do
return
end
subroutine resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
!*****************************************************************************80
!
!! resid() calculates residual contributions to thermodynamic quantities.
!
! Discussion:
!
! The residual function consists of 40 terms. The first 36 are
! used in a global least squares fit to experimental data.
!
! Three terms were added that contribute only in the immediate
! neighborhood of the critical point
! (tk-5) <= T <= (tk+5) C
! 0.20 <= rho <= 0.44 g/cm3,
!
! A single term was added for the region of high pressure and
! low temperature: T < 75 C, 300 MPa < P.
!
! Except in these limited regions, the residual function is
! given by the first 36 terms. The equation is
!
! A(residual)(rho,T)=
! sum(i=1 to 36) (g(i)/k(i)) * (T0/T)**(l(i)) (1-exp(-rho))^(k(i))
! + sum(i=37 to 40) g(i)*delta(i)^(k(i))
! * exp(-alpha(i)*delta(i)^(k(i)) - beta(i)*tau(i)^2)
! (Equation 5)
!
! where
!
! g(i) are coefficients determined by fits to data,
! delta(i) are reduced densities (delta(i)=((rho-rho(i))/rho(i))
! tau(i) are reduced temperatures (tau(i)=((T-tau(i))/tau(i))
! rho(i) are specified densities.
! tau(i) are specified temperatures.
! The k(i) and l(i) are specified integers.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 22 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO, the density, in G/CM3.
!
! Output, real ( kind = rk8 ) AR, the residual contribution to the
! Helmholtz function, in KJ/kg.
!
! Output, real ( kind = rk8 ) CVR, the residual contribution to the
! isochoric (constant volume) heat capacity, in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) DPDRR, the residual contribution to
! the partial derivative dP(T,RHO)/dRHO, with T held fixed, in
! (MegaPascals CM3)/G.
!
! Output, real ( kind = rk8 ) DPDTR, the residual contribution to
! the partial derivative dP(T,RHO)/dT, with RHO held fixed,
! in MegaPascals/degrees Kelvin.
!
! Output, real ( kind = rk8 ) GR, the residual contribution to the Gibbs
! function, in KJ/kg.
!
! Output, real ( kind = rk8 ) HR, the residual contribution to the
! enthalpy, in KJ/kg.
!
! Output, real ( kind = rk8 ) PR, the residual contribution to the pressure,
! in MegaPascals.
!
! Output, real ( kind = rk8 ) SR, the residual contribution to the entropy,
! in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) UR, the residual contribution to the
! internal energy, in KJ/kg.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ), parameter, dimension ( 4 ) :: aad = (/ &
34.0D+00, 40.0D+00, 30.0D+00, 1050.0D+00 /)
real ( kind = rk8 ), parameter, dimension ( 4 ) :: aat = (/ &
20000.0D+00, 20000.0D+00, 40000.0D+00, 25.0D+00 /)
real ( kind = rk8 ), parameter, dimension ( 4 ) :: adz = (/ &
0.319D+00, 0.319D+00, 0.319D+00, 1.55D+00 /)
real ( kind = rk8 ) ar
real ( kind = rk8 ) att
real ( kind = rk8 ), parameter, dimension ( 4 ) :: atz = (/ &
640.0D+00, 640.0D+00, 641.6D+00, 270.0D+00 /)
real ( kind = rk8 ) cvr
real ( kind = rk8 ) dadt
real ( kind = rk8 ) ddz
real ( kind = rk8 ) del
real ( kind = rk8 ) dex
real ( kind = rk8 ) dfdt
real ( kind = rk8 ) dpdrr
real ( kind = rk8 ) dpdtr
real ( kind = rk8 ) e
real ( kind = rk8 ) errtol
real ( kind = rk8 ) ex0
real ( kind = rk8 ) ex1
real ( kind = rk8 ) ex2
real ( kind = rk8 ) fct
real ( kind = rk8 ), parameter, dimension ( 40 ) :: g = (/ &
-530.62968529023D+00, 0.22744901424408D+04, 0.78779333020687D+03, &
-69.830527374994D+00, 0.17863832875422D+05,-0.39514731563338D+05, &
0.33803884280753D+05, -0.13855050202703D+05,-0.25637436613260D+06, &
0.48212575981415D+06, -0.34183016969660D+06, 0.12223156417448D+06, &
0.11797433655832D+07, -0.21734810110373D+07, 0.10829952168620D+07, &
-0.25441998064049D+06, -0.31377774947767D+07, 0.52911910757704D+07, &
-0.13802577177877D+07, -0.25109914369001D+06, 0.46561826115608D+07, &
-0.72752773275387D+07, 0.41774246148294D+06, 0.14016358244614D+07, &
-0.31555231392127D+07, 0.47929666384584D+07, 0.40912664781209D+06, &
-0.13626369388386D+07, 0.69625220862664D+06,-0.10834900096447D+07, &
-0.22722827401688D+06, 0.38365486000660D+06, 0.68833257944332D+04, &
0.21757245522644D+05, -0.26627944829770D+04,-0.70730418082074D+05, &
-0.225D+00, -1.68D+00, 0.055D+00, -93.0D+00 /)
real ( kind = rk8 ) gascon
real ( kind = rk8 ) gr
real ( kind = rk8 ) hr
integer i
integer, parameter, dimension ( 40 ) :: ii = (/ &
0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6, &
8,8,8,8,2,2,0,4,2,2,2,4 /)
integer j
integer, parameter, dimension ( 40 ) :: jj = (/ &
2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,&
2,3,5,7,1,4,4,4,0,2,0,0 /)
integer k
integer l
integer nc
real ( kind = rk8 ) pr
real ( kind = rk8 ) q10
real ( kind = rk8 ) q20
real ( kind = rk8 ) q2a
real ( kind = rk8 ) q5t
real ( kind = rk8 ) qm
real ( kind = rk8 ) qp
real ( kind = rk8 ) qr(11)
real ( kind = rk8 ) qt(10)
real ( kind = rk8 ) rho
real ( kind = rk8 ) sr
real ( kind = rk8 ), parameter :: s_ref = 7.6180720166752D+00
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_ref = 647.073D+00
real ( kind = rk8 ) tau
real ( kind = rk8 ) tx
real ( kind = rk8 ), parameter :: u_ref = - 4328.4549774261D+00
real ( kind = rk8 ) ur
real ( kind = rk8 ) v
errtol = sqrt ( epsilon ( errtol ) )
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'RESID - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'RESID - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was RHO = ', rho
stop
end if
nc = 36
dpdrr = 0.0D+00
pr = 0.0D+00
ar = 0.0D+00
dadt = 0.0D+00
cvr = 0.0D+00
dpdtr = 0.0D+00
ex0 = - rho
ex0 = max ( ex0, -225.0D+00 )
ex0 = min ( ex0, 225.0D+00 )
e = exp ( ex0 )
q10 = rho * rho * e
q20 = 1.0D+00 - e
qr(1) = 0.0D+00
qr(2) = q10
do i = 2, 10
qr(i+1) = qr(i) * q20
end do
v = t_ref / t
qt(1) = t / t_ref
do i = 2, 10
qt(i) = qt(i-1) * v
end do
do i = 1, nc
k = ii(i) + 1
l = jj(i)
qp = g(i) * qr(k+1) * qt(l+1)
pr = pr + qp
dpdrr = dpdrr + ( 2.0D+00 / rho - ( 1.0D+00 - e * dble ( k - 1 ) / &
( 1.0D+00 - e ) ) ) * qp
ar = ar + g(i) * qr(k+2) * qt(l+1) / ( rho**2 * e * dble ( k ) &
* gascon ( ) * t )
dfdt = ( 1.0D+00 - e )**k * dble ( 1 - l ) * qt(l+2) / t_ref / dble ( k )
dadt = dadt + g(i) * dfdt
dpdtr = dpdtr + g(i) * dfdt * rho**2 * e * dble ( k ) / ( 1.0D+00 - e )
cvr = cvr + g(i) * dble ( l ) * dfdt / gascon()
end do
qp = 0.0D+00
q2a = 0.0D+00
do j = 37, 40
k = ii(j)
ddz = adz(j-36)
del = rho / ddz - 1.0D+00
if ( abs ( del ) < errtol ) then
del = errtol
end if
ex1 = - aad(j-36) * del**k
ex1 = max ( ex1, -225.0D+00 )
ex1 = min ( ex1, 225.0D+00 )
dex = exp ( ex1 ) * del**jj(j)
att = aat(j-36)
tx = atz(j-36)
tau = ( t / tx ) - 1.0D+00
ex2 = - att * tau * tau
ex2 = max ( ex2, -225.0D+00 )
ex2 = min ( ex2, 225.0D+00 )
q10 = dex * exp ( ex2 )
qm = dble ( jj(j) ) / del - dble ( k ) * aad(j-36) * del**(k-1)
fct = qm * rho**2 * q10 / ddz
q5t = fct * ( 2.0D+00 / rho + qm / ddz ) - ( rho / ddz )**2 * q10 * &
( dble ( jj(j) ) / del**2 + dble ( k * ( k - 1 ) ) * aad(j-36) * &
del**(k-2) )
dpdrr = dpdrr + q5t * g(j)
qp = qp + g(j) * fct
dadt = dadt - 2.0D+00 * g(j) * att * tau * q10 / tx
dpdtr = dpdtr - 2.0D+00 * g(j) * att * tau * fct / tx
q2a = q2a + t * g(j) * att * ( 4.0D+00 * ex2 + 2.0D+00 ) * q10 / tx**2
ar = ar + q10 * g(j) / ( gascon() * t )
end do
cvr = cvr + q2a / gascon()
pr = pr + qp
sr = - dadt / gascon()
ur = ar + sr
!
! Assign dimensions.
!
ar = gascon() * t * ar
cvr = gascon() * cvr
sr = gascon() * sr
ur = gascon() * t * ur
!
! Adjust energies.
!
ar = ar + gascon ( ) * t * s_ref - gascon ( ) * u_ref
sr = sr - gascon ( ) * s_ref
ur = ur - gascon ( ) * u_ref
gr = ar + pr / rho
hr = ur + pr / rho
return
end
subroutine secvir ( t, vir )
!*****************************************************************************80
!
!! secvir() calculates the second virial coefficient at a given temperature.
!
! Discussion:
!
! The second virial coefficient VIR appears in the first correction term
! to the ideal gas equation of state:
!
! P = R * T / volume + VIR / volume^2 + ...
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 28 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) VIR, the second virial coefficient, in CM3/G.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) b1
real ( kind = rk8 ) b1t
real ( kind = rk8 ) b1tt
real ( kind = rk8 ) b2
real ( kind = rk8 ) b2t
real ( kind = rk8 ) b2tt
real ( kind = rk8 ), parameter, dimension ( 5 ) :: g = (/ &
-0.53062968529023D+03, 0.22744901424408D+04, -0.26627944829770D+04, &
0.78779333020687D+03, -0.69830527374994D+02 /)
real ( kind = rk8 ) gascon
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_ref = 647.073D+00
real ( kind = rk8 ) v
real ( kind = rk8 ) vir
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SECVIR - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
call bb ( t, b1, b2, b1t, b2t, b1tt, b2tt )
v = t_ref / t
vir = b2 + ( &
v * ( g(1) &
+ v * ( g(2) &
+ v * ( g(3) &
+ v * ( g(4) &
+ v * v * g(5) ))))) &
/ ( gascon ( ) * t )
return
end
subroutine secvir_values ( n, tc, vir )
!*****************************************************************************80
!
!! secvir_values() returns some values of the second virial coefficient.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 03 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, pages 24-25.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) VIR, the second virial coefficient, in
! m^3/kg.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 19
integer n
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 5.0D+00, 10.0D+00, 20.0D+00, 30.0D+00, &
40.0D+00, 60.0D+00, 90.0D+00, 120.0D+00, 150.0D+00, &
180.0D+00, 210.0D+00, 240.0D+00, 300.0D+00, 400.0D+00, &
500.0D+00, 700.0D+00, 1000.0D+00, 2000.0D+00 /)
real ( kind = rk8 ) vir
real ( kind = rk8 ), save, dimension ( n_data ) :: virvec = (/ &
-98.96D+00, -90.08D+00, -82.29D+00, -69.36D+00, -59.19D+00, &
-51.07D+00, -39.13D+00, -27.81D+00, -20.83D+00, -16.21D+00, &
-12.98D+00, -10.63D+00, -8.85D+00, -6.39D+00, -4.03D+00, &
-2.71D+00, -1.32D+00, -0.39D+00, 0.53D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
vir = 0.0D+00
else
n = n + 1
tc = tcvec(n)
vir = virvec(n)
end if
return
end
subroutine sound ( t, p, c )
!*****************************************************************************80
!
!! sound() computes the speed of sound given temperature and pressure.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 22 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) P, the pressure, in MegaPascals.
!
! Output, real ( kind = rk8 ) C, the speed of sound, in m/s.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) ab
real ( kind = rk8 ) ai
real ( kind = rk8 ) ar
real ( kind = rk8 ) c
real ( kind = rk8 ) cp
real ( kind = rk8 ) cpi
real ( kind = rk8 ) cv
real ( kind = rk8 ) cvb
real ( kind = rk8 ) cvi
real ( kind = rk8 ) cvr
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdrb
real ( kind = rk8 ) dpdrr
real ( kind = rk8 ) dpdt
real ( kind = rk8 ) dpdtb
real ( kind = rk8 ) dpdtr
real ( kind = rk8 ) gb
real ( kind = rk8 ) gi
real ( kind = rk8 ) gr
real ( kind = rk8 ) hb
real ( kind = rk8 ) hi
real ( kind = rk8 ) hr
real ( kind = rk8 ) p
real ( kind = rk8 ) pb
real ( kind = rk8 ) pr
real ( kind = rk8 ) rho
real ( kind = rk8 ) rhol
real ( kind = rk8 ) rho_start
real ( kind = rk8 ) rhov
real ( kind = rk8 ) sb
real ( kind = rk8 ) si
real ( kind = rk8 ) sr
real ( kind = rk8 ) t
real ( kind = rk8 ) t2
real ( kind = rk8 ) ub
real ( kind = rk8 ) ui
real ( kind = rk8 ) ur
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SOUND - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative pressures.
!
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SOUND - Fatal error!'
write ( *, '(a)' ) ' The input pressure P must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was P = ', p
stop
end if
!
! For the given pressure, compute the saturation temperature.
!
rhol = 0.0D+00
rhov = 0.0D+00
call tsat ( p, t2, rhol, rhov )
!
! Depending on whether the temperature is above or below the
! saturation temperature, we expect to compute the density of
! a liquid or vapor.
!
if ( t < t2 ) then
rho_start = 1.9D+00
else
rho_start = 0.01D+00
end if
call dense ( p, t, rho_start, rho, dpdr )
!
! From T and RHO, compute the thermodynamic properties.
!
call ideal ( t, ai, cpi, cvi, gi, hi, si, ui )
call resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
call base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
cv = cvb + cvr + cvi
dpdr = dpdrb + dpdrr
dpdt = dpdtb + dpdtr
cp = cv + t * dpdt * dpdt / ( dpdr * rho**2 )
c = sqrt ( 1000.0D+00 * cp * dpdr / cv )
return
end
subroutine sound_values ( n, tc, p, c )
!*****************************************************************************80
!
!! sound_values() returns some values of the speed of sound.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, page 238-246.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) C, the speed of sound, in m/s.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 20
real ( kind = rk8 ) c
real ( kind = rk8 ), save, dimension ( n_data ) :: cvec = (/ &
1401.0D+00, 472.8D+00, 533.7D+00, 585.7D+00, 609.5D+00, &
632.2D+00, 674.6D+00, 713.9D+00, 802.0D+00, 880.1D+00, &
1017.8D+00, 1115.9D+00, 1401.7D+00,1402.6D+00, 1409.6D+00, &
1418.1D+00, 1443.1D+00, 1484.6D+00, 1577.1D+00, 1913.4D+00 /)
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 5.0D+00, 10.0D+00, 50.0D+00, &
100.0D+00, 250.0D+00, 500.0D+00, 1000.0D+00, 2500.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 100.0D+00, 200.0D+00, 300.0D+00, 350.0D+00, &
400.0D+00, 500.0D+00, 600.0D+00, 850.0D+00, 1100.0D+00, &
1600.0D+00, 2000.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
c = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
c = cvec(n)
end if
return
end
subroutine surten ( t, sigma )
!*****************************************************************************80
!
!! surten() returns the surface tension as a function of temperature.
!
! Discussion:
!
! SURTEN uses an equation that yields values of the surface tension to
! within the accuracy of measurements from the triple point to the
! critical point.
!
! Sigma = B * ( (TSTAR-T)/TSTAR)^Mu * (1+b*(TSTAR-T)/TSTAR)
!
! where:
!
! TSTAR = 647.15 Degrees Kelvin,
! B = 0.2358 Pascals * Meters
! b = -0.625,
! Mu = 1.256.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 13 January 2013
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) SIGMA, the surface tension,
! in Pascal * m = Newton / m.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ), parameter :: b_cap = 0.2358D+00
real ( kind = rk8 ), parameter :: b_small = -0.625D+00
real ( kind = rk8 ), parameter :: mu = 1.256D+00
real ( kind = rk8 ) sigma
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_star = 647.15D+00
real ( kind = rk8 ) term
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SURTEN - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
term = ( t_star - t ) / t_star
sigma = b_cap * term**mu * ( 1.0D+00 + b_small * term )
!
! Need this conversion to match the table, but justification is there none.
!
sigma = 1000.0D+00 * sigma
return
end
subroutine surten_values ( n, tc, sigma )
!*****************************************************************************80
!
!! surten_values() returns some values of the surface tension.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, pages 267.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) SIGMA, the surface tension,
! in Pascal * m = Newton / m.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 14
integer n
real ( kind = rk8 ) sigma
real ( kind = rk8 ), save, dimension ( n_data ) :: sigmavec = (/ &
74.22D+00, 72.74D+00, 71.20D+00, 69.60D+00, 67.95D+00, &
58.92D+00, 48.75D+00, 37.68D+00, 26.05D+00, 14.37D+00, &
8.78D+00, 3.67D+00, 0.40D+00, 0.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
10.0D+00, 20.0D+00, 30.0D+00, 40.0D+00, 50.0D+00, &
100.0D+00, 150.0D+00, 200.0D+00, 250.0D+00, 300.0D+00, &
325.0D+00, 350.0D+00, 370.0D+00, 373.976D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
sigma = 0.0D+00
else
n = n + 1
tc = tcvec(n)
sigma = sigmavec(n)
end if
return
end
subroutine tdpsdt ( t, dp )
!*****************************************************************************80
!
!! tdpsdt() computes the quantity T * dP(Sat)/dT.
!
! Discussion:
!
! Here T is the temperature and P(Sat) is the vapor pressure.
! It is used by TSAT_EST and TSAT.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 13 January 2013
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) DP, the value T*(dP(Sat)/dT),
! in MegaPascals.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ), parameter, dimension ( 8 ) :: a = (/ &
-7.8889166D+00, 2.5514255D+00, -6.716169D+00, 33.239495D+00, &
-105.38479D+00, 174.35319D+00, -148.39348D+00, 48.631602D+00 /)
real ( kind = rk8 ) b
real ( kind = rk8 ) c
real ( kind = rk8 ) dp
integer i
real ( kind = rk8 ) q
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_ref = 647.25D+00
real ( kind = rk8 ) v
real ( kind = rk8 ) w
real ( kind = rk8 ) y
real ( kind = rk8 ) z
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TDPSDT - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
v = t / t_ref
w = 1.0D+00 - v
b = 0.0D+00
c = 0.0D+00
do i = 1, 8
z = dble ( i + 1 ) / 2.0D+00
y = a(i) * w**z
c = c + ( y / w ) * ( 0.5D+00 - 0.5D+00 * dble ( i ) - 1.0D+00 / v )
b = b + y
end do
q = b / v
dp = 22.093D+00 * exp ( q ) * c
return
end
subroutine thercon ( t, rho, lambda )
!*****************************************************************************80
!
!! thercon() calculates thermal conductivity for given temperature and density.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 20 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO, the density, in G/CM3.
!
! Output, real ( kind = rk8 ) LAMBDA, the thermal conductivity,
! in mW/(m degrees Kelvin).
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) a
real ( kind = rk8 ), parameter, dimension ( 0:3 ) :: acof = (/ &
2.02223D+00, 14.11166D+00, 5.25597D+00, -2.01870D+00 /)
real ( kind = rk8 ), parameter :: a_con = 18.66D+00
real ( kind = rk8 ) b(0:4,0:5)
real ( kind = rk8 ), parameter :: b_con = 1.00D+00
real ( kind = rk8 ), parameter :: c_con = 3.7711D-08
real ( kind = rk8 ) chi
real ( kind = rk8 ) cjth
real ( kind = rk8 ) cjtt
real ( kind = rk8 ) cp
real ( kind = rk8 ) cv
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdr2
real ( kind = rk8 ) dpdt
real ( kind = rk8 ) eta
real ( kind = rk8 ) g
real ( kind = rk8 ) h
integer i
integer j
real ( kind = rk8 ) lambda
real ( kind = rk8 ) lambda0
real ( kind = rk8 ) lambda_del
real ( kind = rk8 ), parameter :: omega = 0.4678D+00
real ( kind = rk8 ) p
real ( kind = rk8 ), parameter :: p_ref = 22.115D+00
real ( kind = rk8 ) power
real ( kind = rk8 ) rho
real ( kind = rk8 ), parameter :: rho_ref = 317.763D+00
real ( kind = rk8 ) rho2
real ( kind = rk8 ) s
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_ref = 647.27D+00
real ( kind = rk8 ) total
real ( kind = rk8 ) u
data b / &
1.3293046D+00, 1.7018363D+00, 5.2246158D+00, &
8.7127675D+00, -1.8525999D+00, &
-0.40452437D+00, -2.2156845D+00, -10.124111D+00, &
-9.5000611D+00, 0.93404690D+00, &
0.24409490D+00, 1.6511057D+00, 4.9874687D+00, &
4.3786606D+00, 0.0D+00, &
0.018660751D+00, -0.76736002D+00, -0.27297694D+00, &
-0.91783782D+00, 0.0D+00, &
-0.12961068D+00, 0.37283344D+00, -0.43083393D+00, &
0.0D+00, 0.0D+00, &
0.044809953D+00, -0.11203160D+00, 0.13333849D+00, &
0.0D+00, 0.0D+00 /
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'THERCON - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'THERCON - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was RHO = ', rho
stop
end if
!
! Compute DPDR, DPDT, ETA.
!
call therm ( t, rho, a, cjth, cjtt, cp, cv, dpdr, dpdt, g, h, p, s, u )
call viscosity ( t, rho, eta )
!
! Convert RHO from G/CM3 to kg/M3,
! Convert DPDR from ? to ?.
!
rho2 = 1000.0D+00 * rho
dpdr2 = dpdr / 1000.0D+00
!
! Compute LAMBDA0.
!
total = 0.0D+00
do i = 0, 3
total = total + acof(i) * ( t_ref / t )**i
end do
lambda0 = sqrt ( t / t_ref ) / total
!
! Compute CHI.
!
chi = rho2 * p_ref / ( rho_ref**2 * dpdr2 )
!
! Compute delta_Lambda.
!
power = - a_con * ( ( t_ref - t ) / t )**2 - b_con * ( ( rho2 - rho_ref ) &
/ rho_ref )**4
lambda_del = ( c_con / eta ) * ( ( t * rho_ref ) / ( t_ref * rho ) )**2 &
* ( t_ref / p_ref )**2 * dpdt**2 * chi**omega * sqrt ( rho2 / rho_ref ) &
* exp ( power )
!
! Compute LAMBDA.
!
total = 0.0D+00
do i = 0, 4
do j = 0, 5
total = total + b(i,j) * ( ( t_ref - t ) / t )**i * &
( ( rho2 - rho_ref ) / rho_ref )**j
end do
end do
lambda = lambda0 * exp ( ( rho2 / rho_ref ) * total ) + lambda_del
!
! Temporary fix.
!
lambda = 1000.0D+00 * lambda
return
end
subroutine thercon_values ( n, tc, p, lambda )
!*****************************************************************************80
!
!! thercon_values() returns some values of the thermal conductivity.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, page 264.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) LAMBDA, the thermal conductivity, in
! mW/(m degrees Kelvin).
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 35
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ) lambda
real ( kind = rk8 ), save, dimension ( n_data ) :: lambdavec = (/ &
561.0D+00, 561.3D+00, 561.5D+00, 562.4D+00, 563.7D+00, &
565.1D+00, 566.5D+00, 567.9D+00, 569.3D+00, 570.6D+00, &
572.0D+00, 573.4D+00, 574.8D+00, 576.1D+00, 577.5D+00, &
580.2D+00, 582.9D+00, 585.5D+00, 588.1D+00, 590.7D+00, &
593.3D+00, 595.8D+00, 598.3D+00, 603.1D+00, 607.8D+00, &
612.2D+00, 607.2D+00, 643.6D+00, 666.8D+00, 25.08D+00, &
28.85D+00, 33.28D+00, 54.76D+00, 79.89D+00, 107.3D+00 /)
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
1.0D+00, 5.0D+00, 10.0D+00, 25.0D+00, 50.0D+00, &
75.0D+00, 100.0D+00, 125.0D+00, 150.0D+00, 175.0D+00, &
200.0D+00, 225.0D+00, 250.0D+00, 275.0D+00, 300.0D+00, &
350.0D+00, 400.0D+00, 450.0D+00, 500.0D+00, 550.0D+00, &
600.0D+00, 650.0D+00, 700.0D+00, 800.0D+00, 900.0D+00, &
1000.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 25.0D+00, 50.0D+00, 75.0D+00, 100.0D+00, &
150.0D+00, 200.0D+00, 400.0D+00, 600.0D+00, 800.0D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
lambda = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
lambda = lambdavec(n)
end if
return
end
subroutine therm ( t, rho, a, cjth, cjtt, cp, cv, dpdr, dpdt, g, h, p, s, u )
!*****************************************************************************80
!
!! therm() calculates thermodynamic functions given temperature and density.
!
! Discussion:
!
! Thermodynamic values were calculated from an analytic equation
! that approximates the Helmholtz function (specific Helmholtz
! energy) for ordinary water and steam, of the form A=A(rho,T)
! where A is the Helmholtz function, rho the density, and T
! the absolute (thermodynamic) temperature. Any thermodynamic
! value for any state, liquid, vapor or metastable, may be
! calculated by differentiation of this equation in accord with
! the first and second laws of thermodynamics.
!
! The International Association for the Properties of Steam
! has provisionally accepted this formulation for the range
! 273.15 <= T <= 1273.15 degrees Kelvin, where, for 423.15 <= T,
! the maximum pressure is Pmax = 1500 MPa = 15000 bar, and for
! 273.15 <= T < 423.15, the maximum pressure is
! Pmax = 100 * (5 + (T-273.15)/15) MPa.
!
! Close to the critical point, a small region is excluded:
! Abs(T-Tk) < 1, abs((rho-rhok)/rhok) < 0.3.
!
! The equation has a wider useful range, namely, fluid states
! of pure, undissociated water and steam defined by
! 260 <= T <= 2500 K and 0 <= P <= 3000 MPa.
!
! Thermodynamic property values for specific volume, density,
! specific internal energy, specific enthalpy, and specific
! entropy of water and steam were tabulated over the range
! 0 <= t <= 2000 C, 0 <= P <= 3000 MPa. The reference
! state is the liquid at the triple point, for which the
! internal energy and entropy have been assigned the value zero.
!
! Thermodynamic quantities are determined from the Helmholtz function
! A(rho,T), which is computed as the sum of three terms:
!
! A(rho,T) = A(base)(rho,T) + A(residual)(rho,T) + A(ideal)(T)
! (Equation 1)
!
! Because A(rho,T) is everywhere single valued and analytic,
! we can derive closed form relations for all other properties.
! In the following, unless otherwise indicated, the independent
! variables are temperature T and density RHO, and differentiation
! with respect to one variable is to imply that the other is fixed.
!
! Pressure: P = RHO^2 * dA/dRHO
! Density derivative: dP/dRHO = 2*P/RHO + RHO^2 * d2A/dRHO2
! Temperature derivative: dP/dT = RHO^2 * d2A/(dRHO dT)
! Specific entropy: S = - dA/dT
! Specific internal energy: U = A + T*S
! Specific enthalpy: H = U + P/RHO
! Specific heat capacity
! at constant volume: Cv = - T * d2A/dT2
! Specific Gibbs function: G = A + P/RHO
! Specific heat capacity
! at constant pressure: Cp = Cv + (T*(dP/dT)^2)/(RHO^2*dP/dRHO)
! Speed of sound: Omega = Sqrt ((Cp/Cv) * dP/dRHO)
! Second virial coefficient: B = 1/(2*R*T) * (d2P/dRHO2) (at RHO=0)
! Isothermal Joule-Thomson
! coefficient: DeltaT = (dH/dP) (fixed T) =
! (1/RHO)-(T*dP/dT)/(RHO^2*dP/dRHO)
! Joule-Thomson coefficient: Mu = (dT/dP) (fixed H) = DeltaT/Cp
! Isentropic temperature-
! pressure coefficient: BetaS = (dT/dP) (fixed S) =
! (DeltaT - 1/RHO)/Cp
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 19 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO, the fluid density, in G/CM3.
!
! Output, real ( kind = rk8 ) A, the Helmholtz function, in KJ/kg.
!
! Output, real ( kind = rk8 ) CJTH, the Joule-Thomson coefficient,
! in K/MegaPascals.
!
! Output, real ( kind = rk8 ) CJTT, the isothermal Joule-Thomson coefficient,
! in CM3/G.
!
! Output, real ( kind = rk8 ) CP, the isobaric (constant pressure) heat
! capacity, in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) CV, the isochoric (constant volume) heat
! capacity, in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) DPDR, the partial derivative
! dP(T,RHO)/dRHO, with T held fixed, in MegaPascals*CM3/G.
!
! Output, real ( kind = rk8 ) DPDT, the partial derivative
! dP(T,RHO)/dT, with RHO held fixed, in MegaPascals/degrees Kelvin.
!
! Output, real ( kind = rk8 ) G, the Gibbs free energy, in KJ/kg.
!
! Output, real ( kind = rk8 ) H, the enthalpy, in KJ/kg.
!
! Output, real ( kind = rk8 ) P, the pressure, in MegaPascals.
!
! Output, real ( kind = rk8 ) S, the entropy, in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) U, the internal energy, in KJ/kg.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) a
real ( kind = rk8 ) ab
real ( kind = rk8 ) ai
real ( kind = rk8 ) ar
real ( kind = rk8 ) cjth
real ( kind = rk8 ) cjtt
real ( kind = rk8 ) cp
real ( kind = rk8 ) cpi
real ( kind = rk8 ) cv
real ( kind = rk8 ) cvb
real ( kind = rk8 ) cvi
real ( kind = rk8 ) cvr
logical, parameter :: debug = .false.
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdrb
real ( kind = rk8 ) dpdrr
real ( kind = rk8 ) dpdt
real ( kind = rk8 ) dpdtb
real ( kind = rk8 ) dpdtr
real ( kind = rk8 ) g
real ( kind = rk8 ) gb
real ( kind = rk8 ) gi
real ( kind = rk8 ) gr
real ( kind = rk8 ) h
real ( kind = rk8 ) hb
real ( kind = rk8 ) hi
real ( kind = rk8 ) hr
real ( kind = rk8 ) p
real ( kind = rk8 ) pb
real ( kind = rk8 ) pr
real ( kind = rk8 ) rho
real ( kind = rk8 ) s
real ( kind = rk8 ) sb
real ( kind = rk8 ) si
real ( kind = rk8 ) sr
real ( kind = rk8 ) t
real ( kind = rk8 ) u
real ( kind = rk8 ) ub
real ( kind = rk8 ) ui
real ( kind = rk8 ) ur
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'THERM - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'THERM - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was RHO = ', rho
stop
end if
call ideal ( t, ai, cpi, cvi, gi, hi, si, ui )
call resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
call base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
a = ab + ar + ai
cv = cvb + cvr + cvi
if ( debug ) then
write ( *, * ) ' '
write ( *, * ) 'THERM:'
write ( *, * ) ' CVB = ', cvb
write ( *, * ) ' CVR = ', cvr
write ( *, * ) ' CVI = ', cvi
write ( *, * ) ' CV = ', cv
end if
dpdr = dpdrb + dpdrr
dpdt = dpdtb + dpdtr
p = pb + pr
s = sb + sr + si
u = ub + ur + ui
if ( debug ) then
write ( *, * ) ' '
write ( *, * ) 'THERM:'
write ( *, * ) ' UB = ', ub
write ( *, * ) ' UR = ', ur
write ( *, * ) ' UI = ', ui
end if
g = a + p / rho
h = u + p / rho
cp = cv + t * dpdt**2 / ( dpdr * rho**2 )
cjtt = 1.0D+00 / rho - t * dpdt / ( dpdr * rho**2 )
cjth = - cjtt / cp
return
end
subroutine tsat ( p, t, rhol, rhov )
!*****************************************************************************80
!
!! tsat() calculates the saturation temperature for a given pressure.
!
! Discussion:
!
! The corresponding liquid and vapor densities are also computed.
! The saturation temperature is also known as the "vapor temperature".
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) P, the vapor pressure, in MegaPascals.
!
! Output, real ( kind = rk8 ) T, the vapor temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) RHOL, the liquid density, in G/CM3.
!
! Output, real ( kind = rk8 ) RHOV, the vapor density, in G/CM3.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
logical, parameter :: debug = .false.
real ( kind = rk8 ) delg
real ( kind = rk8 ) dp
real ( kind = rk8 ) dp2
real ( kind = rk8 ) errtol
real ( kind = rk8 ) gascon
integer it
integer, parameter :: it_max = 50
real ( kind = rk8 ) p
real ( kind = rk8 ) p_consistent
real ( kind = rk8 ) rhol
real ( kind = rk8 ) rhov
real ( kind = rk8 ) t
!
! Initialize output quantities, obliterating any input value.
!
t = 0.0D+00
rhol = 0.0D+00
rhov = 0.0D+00
!
! Set the error tolerance.
!
errtol = sqrt ( epsilon ( errtol ) )
!
! Refuse to handle zero or negative pressure.
!
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TSAT - Fatal error!'
write ( *, '(a)' ) ' The input pressure must be positive!'
write ( *, '(a,g14.6)' ) ' Your value was P = ', p
stop
end if
!
! Estimate the saturation temperature.
!
call tsat_est ( p, t )
if ( debug ) then
write ( *, * ) ' '
write ( *, * ) 'TSAT:'
write ( *, '(2g14.6)' ) p, t, rhol, rhov
end if
do it = 1, it_max
call corr ( t, p, p_consistent, rhol, rhov, delg )
dp = delg * gascon ( ) * t * rhol * rhov / ( rhol - rhov )
call tdpsdt ( t, dp2 )
t = t * ( 1.0D+00 - dp / dp2 )
if ( debug ) then
write ( *, '(2g14.6)' ) p, t, rhol, rhov
end if
if ( abs ( delg ) < errtol ) then
return
end if
end do
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TSAT - Warning!'
write ( *, '(a)' ) ' The iteration did not converge.'
write ( *, '(a,i6)' ) ' Number of iterations was ', it_max
write ( *, '(a,g14.6)' ) ' Last iterate was ', t
write ( *, '(a,g14.6)' ) ' Last DELG was ', delg
return
end
subroutine tsat_est ( p, t )
!*****************************************************************************80
!
!! tsat_est() makes a rough estimate of the saturation temperature.
!
! Discussion:
!
! The saturation temperature is also called the vapor temperature.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 February 2002
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) P, the pressure, in MegaPascals. The tabulated
! range for P is
! 0.00061173 MegaPascals <= P <= P_CRIT = 22.055 MegaPascals.
! The input value of P must be positive.
!
! Output, real ( kind = rk8 ) T, the saturation temperature,
! in degrees Kelvin. This value will always be in the range
! [ 273.15, 647.126 ].
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: npol = 4
real ( kind = rk8 ), parameter, dimension ( 0:npol ) :: c = (/ &
372.83D+00, 27.7589D+00, 2.3819D+00, 0.24834D+00, 0.0193855D+00 /)
real ( kind = rk8 ) dp
real ( kind = rk8 ) dt
real ( kind = rk8 ) errtol
integer it
integer, parameter :: it_max = 8
real ( kind = rk8 ) p
real ( kind = rk8 ) pl
real ( kind = rk8 ), parameter :: p_crit = 22.055D+00
real ( kind = rk8 ) pp
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_crit = 647.126D+00
real ( kind = rk8 ), parameter :: t_min = 273.15D+00
real ( kind = rk8 ) t_old
errtol = sqrt ( epsilon ( errtol ) )
!
! Refuse to handle zero or negative pressure.
!
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TSAT_EST - Fatal error!'
write ( *, '(a)' ) ' The input pressure must be positive!'
write ( *, '(a,g14.6)' ) ' Your value was P = ', p
stop
end if
if ( p_crit < p ) then
t = t_crit
return
end if
!
! The initial estimate for T uses a polyonmial in the logarithm of P.
!
pl = 2.302585D+00 + log ( p )
call r8poly_val_horner ( npol, c, pl, t )
t = min ( t, t_crit )
t = max ( t, t_min )
dt = 0.0D+00
do it = 1, it_max
call psat_est ( t, pp )
call tdpsdt ( t, dp )
if ( abs ( p - pp ) < errtol * p ) then
return
end if
dt = t * ( p - pp ) / dp
t_old = t
t = t * ( 1.0D+00 + ( p - pp ) / dp )
t = min ( t, t_crit )
t = max ( t, t_min )
if ( abs ( dt ) < errtol * ( abs ( t ) + 1.0D+00 ) ) then
return
else if ( abs ( t - t_old ) < errtol ) then
return
end if
end do
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TSAT_EST - Warning!'
write ( *, '(a)' ) ' The iteration did not converge.'
write ( *, '(a,i6)' ) ' Number of iterations was ', it_max
write ( *, '(a,g14.6)' ) ' Convergence tolerance was ', errtol
write ( *, '(a,g14.6)' ) ' Last iterate was ', t
write ( *, '(a,g14.6)' ) ' Last correction was ', dt
return
end
subroutine tsat_values ( n, p, tc )
!*****************************************************************************80
!
!! tsat_values() returns some values of the saturation temperature.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 05 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, pages 16-22.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) TC, the saturation temperature, in
! degrees Celsius.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 20
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
0.0061173D+00, 0.012D+00, 0.025D+00, 0.055D+00, 0.080D+00, &
0.11D+00, 0.16D+00, 0.25D+00, 0.50D+00, 0.75D+00, &
1.0D+00, 1.5D+00, 2.0D+00, 5.0D+00, 10.0D+00, &
20.0D+00, 50.0D+00, 100.0D+00, 200.0D+00, 220.55D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.010D+00, 9.655D+00, 21.080D+00, 34.589D+00, 41.518D+00, &
47.695D+00, 55.327D+00, 64.980D+00, 81.339D+00, 91.783D+00, &
99.632D+00, 111.378D+00, 120.443D+00, 151.866D+00, 179.916D+00, &
212.417D+00, 263.977D+00, 311.031D+00, 365.800D+00, 373.976D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
p = 0.0D+00
tc = 0.0D+00
else
n = n + 1
p = pvec(n)
tc = tcvec(n)
end if
return
end
subroutine viscosity ( t, rho, eta )
!*****************************************************************************80
!
!! viscosity() calculates the viscosity for given temperature and density.
!
! Discussion:
!
! On 02 February 2002, I discovered that the Haar/Gallagher/Kell
! reference apparently reversed the sign on the A3 coefficient.
! That made the results better, but still off.
!
! Apparently Haar/Gallagher/Kell had a transcription error in
! the value of B(4,1), which they list as -0.273093, but which
! should be -0.253093.
!
! These two corrections courtesy of Meyer/McClintock/Silvestri/Spencer.
!
! Now the results look proper! And just 12 years late...
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 February 2002
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! International Association for the Properties of Steam,
! Release on Dynamic Viscosity of Water Substance,
! National Bureau of Standards, Washington DC, 1975, revised 1983.
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO, the density, in G/CM3.
!
! Output, real ( kind = rk8 ) ETA, the viscosity, in MegaPascal seconds.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: npol_t = 3
real ( kind = rk8 ), parameter, dimension ( 0:npol_t ) :: a = (/ &
0.0181583D+00, 0.0177624D+00, 0.0105287D+00, -0.0036744D+00 /)
real ( kind = rk8 ) arg
real ( kind = rk8 ), dimension(0:5,0:4) :: b = reshape ( (/ &
0.501938D+00, 0.162888D+00, -0.130356D+00, &
0.907919D+00, -0.551119D+00, 0.146543D+00, &
0.235622D+00, 0.789393D+00, 0.673665D+00, &
1.207552D+00, 0.0670665D+00, -0.0843370D+00, &
-0.274637D+00, -0.743539D+00, -0.959456D+00, &
-0.687343D+00, -0.497089D+00, 0.195286D+00, &
0.145831D+00, 0.263129D+00, 0.347247D+00, &
0.213486D+00, 0.100754D+00, -0.032932D+00, &
-0.0270448D+00, -0.0253093D+00, -0.0267758D+00, &
-0.0822904D+00, 0.0602253D+00, -0.0202595D+00 /), &
(/ 6,5 /) )
logical, parameter :: debug = .false.
real ( kind = rk8 ) eta
real ( kind = rk8 ) eta0
integer i
integer j
real ( kind = rk8 ) rho
real ( kind = rk8 ), parameter :: rho_max = 1.050D+00
real ( kind = rk8 ), parameter :: rho_ref = 0.317763D+00
real ( kind = rk8 ) t
real ( kind = rk8 ) total
real ( kind = rk8 ), parameter :: t_max = 800.00D+00
real ( kind = rk8 ), parameter :: t_ref = 647.27D+00
logical, save :: warning1 = .false.
logical, save :: warning2 = .false.
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'VISCOSITY - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was T = ', t
stop
end if
if ( t_max < t .and. ( .not. warning1 ) ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'VISCOSITY - Warning (once only)!'
write ( *, '(a,g14.6)' ) &
' The input temperature T should be no more than ', t_max
write ( *, '(a,g14.6)' ) ' The input value was T = ', t
write ( *, '(a)' ) ' '
warning1 = .true.
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'VISCOSITY - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was RHO = ', rho
stop
end if
if ( rho_max < rho .and. ( .not. warning2 ) ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'VISCOSITY - Warning (once only)!'
write ( *, '(a,g14.6)' ) &
' The input density RHO should be no more than ', rho_max
write ( *, '(a,g14.6)' ) ' The input value was RHO = ', rho
write ( *, '(a)' ) ' '
warning2 = .true.
end if
!
! Compute ETA0.
!
arg = t_ref / t
call r8poly_val_horner ( npol_t, a, arg, total )
eta0 = sqrt ( t / t_ref ) / total
!
! Compute ETA.
!
total = 0.0D+00
do i = 0, 5
do j = 0, 4
total = total + b(i,j) * ( ( t_ref - t ) / t )**i * &
( ( rho - rho_ref ) / rho_ref )**j
end do
end do
eta = eta0 * exp ( ( rho / rho_ref ) * total )
return
end
subroutine viscosity_values ( n, tc, p, eta )
!*****************************************************************************80
!
!! viscosity_values() returns some values of the viscosity function for testing.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, page 263.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) ETA, the viscosity, in MegaPascal seconds.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 34
real ( kind = rk8 ) eta
real ( kind = rk8 ), save, dimension ( n_data ) :: etavec = (/ &
1792.0D+00, 1791.0D+00, 1790.0D+00, 1786.0D+00, 1780.0D+00, &
1775.0D+00, 1769.0D+00, 1764.0D+00, 1759.0D+00, 1754.0D+00, &
1749.0D+00, 1744.0D+00, 1739.0D+00, 1735.0D+00, 1731.0D+00, &
1722.0D+00, 1714.0D+00, 1707.0D+00, 1700.0D+00, 1694.0D+00, &
1687.0D+00, 1682.0D+00, 1676.0D+00, 1667.0D+00, 1659.0D+00, &
1653.0D+00, 890.8D+00, 547.1D+00, 378.4D+00, 12.28D+00, &
16.18D+00, 24.45D+00, 32.61D+00, 40.38D+00 /)
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
1.0D+00, 5.0D+00, 10.0D+00, 25.0D+00, 50.0D+00, &
75.0D+00, 100.0D+00, 125.0D+00, 150.0D+00, 175.0D+00, &
200.0D+00, 225.0D+00, 250.0D+00, 275.0D+00, 300.0D+00, &
350.0D+00, 400.0D+00, 450.0D+00, 500.0D+00, 550.0D+00, &
600.0D+00, 650.0D+00, 700.0D+00, 800.0D+00, 900.0D+00, &
1000.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 25.0D+00, 50.0D+00, 75.0D+00, 100.0D+00, &
200.0D+00, 400.0D+00, 600.0D+00, 800.0D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
eta = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
eta = etavec(n)
end if
return
end
subroutine volume ( t, rho, v, dvdt, dvdr )
!*****************************************************************************80
!
!! volume() computes specific volume derivatives given temperature and density.
!
! Discussion:
!
! Because A(rho,T) is everywhere single valued and analytic,
! we can derive closed form relations for all other properties.
!
! The independent variables are temperature T and density RHO,
! and differentiation with respect to one variable is to imply that
! the other is held at a fixed value.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 28 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO, the fluid density, in G/CM3.
!
! Output, real ( kind = rk8 ) V, the specific volume, in CM3/G.
!
! Output, real ( kind = rk8 ) DVDT, the partial derivative dV(T,RHO)/dT,
! where V is the specific volume, in CM3 / (G * degrees Kelvin).
!
! Output, real ( kind = rk8 ) DVDR, the partial derivative dV(T,RHO)/dRHO,
! where V is the specific volume, in CM3^2 / ( G^2 ).
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) ab
real ( kind = rk8 ) ar
real ( kind = rk8 ) cvb
real ( kind = rk8 ) cvr
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdrb
real ( kind = rk8 ) dpdrr
real ( kind = rk8 ) dpdt
real ( kind = rk8 ) dpdtb
real ( kind = rk8 ) dpdtr
real ( kind = rk8 ) dvdr
real ( kind = rk8 ) dvdt
real ( kind = rk8 ) gb
real ( kind = rk8 ) gr
real ( kind = rk8 ) hb
real ( kind = rk8 ) hr
real ( kind = rk8 ) pb
real ( kind = rk8 ) pr
real ( kind = rk8 ) rho
real ( kind = rk8 ) sb
real ( kind = rk8 ) sr
real ( kind = rk8 ) t
real ( kind = rk8 ) ub
real ( kind = rk8 ) ur
real ( kind = rk8 ) v
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'VOLUME - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'VOLUME - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was RHO = ', rho
stop
end if
call resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
call base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
dpdr = dpdrb + dpdrr
dpdt = dpdtb + dpdtr
dvdt = dpdt / ( dpdr * rho**2 )
!
! Because V = 1/Rho, dV/dRho = -1/Rho**2
!
dvdr = - 1.0D+00 / rho**2
v = 1.0D+00 / rho
return
end
! =========================================================================
! Source File: steam_calc.f90
! =========================================================================
program steam_calc
implicit none
integer, parameter :: rk8 = kind(1.0D+00)
! Input variables
integer :: input_type
double precision :: val1, val2
! Conversions
double precision :: t_c, t_k, p_kpa, p_mpa
double precision :: x_qual, spec_vol, spec_u, spec_h, spec_s, density
double precision :: cp_val, cv_val, speed_sound, visc_eta, therm_lambda
! Saturated values
double precision :: t_sat_k, p_sat_mpa, rhol, rhov
double precision :: h_liq, h_vap, s_liq, s_vap, u_liq, u_vap
double precision :: cp_liq, cp_vap, cv_liq, cv_vap, sound_liq, sound_vap
double precision :: visc_liq, visc_vap, therm_liq, therm_vap
! Calculation values
character(len=30) :: state_name
double precision :: rho, rho_start, dpdr, dpdt
double precision :: ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub
double precision :: ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur
double precision :: a_term, cjth, cjtt, cp_term, cv_term, g_term, h_term, p_term, s_term, u_term
! Solver variables
double precision :: t_min, t_max, t_mid, h_mid, s_mid, err
integer :: iter
logical :: is_mix
! Functions from steam.f90
double precision :: gascon
! Read from stdin
read(*,*) input_type
read(*,*) val1
read(*,*) val2
! Default quality set to -1 (meaning not in mixture region, single phase)
x_qual = -1.0d0
is_mix = .false.
select case(input_type)
! -------------------------------------------------------------
! 1. Temperature and Pressure (val1 = T [C], val2 = P [kPa])
! -------------------------------------------------------------
case(1)
t_c = val1
t_k = t_c + 273.15d0
p_kpa = val2
p_mpa = p_kpa / 1000.0d0
! Check boundaries
if (t_k < 273.16d0 .or. t_k > 1273.15d0) then
write(*,*) 'ERROR: Temperature must be between 0.01 C and 1000 C'
stop
end if
if (p_mpa < 0.000611d0 .or. p_mpa > 100.0d0) then
write(*,*) 'ERROR: Pressure must be between 0.61 kPa and 100000 kPa'
stop
end if
if (p_mpa < 22.055d0) then
call tsat(p_mpa, t_sat_k, rhol, rhov)
if (abs(t_k - t_sat_k) < 1.0d-3) then
! On the saturation boundary, assume superheated or mixture
! But T & P is not enough to define mixture state, so default to saturated liquid
state_name = 'Saturated boundary'
rho_start = rhol
else if (t_k < t_sat_k) then
state_name = 'Subcooled Liquid'
rho_start = 1.11d0 - 0.0004d0 * t_k
else
state_name = 'Superheated Vapor'
rho_start = p_mpa / (gascon() * t_k)
end if
else
state_name = 'Supercritical Fluid'
if (t_k < 647.13d0) then
rho_start = 1.11d0 - 0.0004d0 * t_k
else
rho_start = p_mpa / (gascon() * t_k)
end if
end if
call dense(p_mpa, t_k, rho_start, rho, dpdr)
call therm(t_k, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, spec_h, p_term, spec_s, spec_u)
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
! -------------------------------------------------------------
! 2. Temperature and Quality (val1 = T [C], val2 = x)
! -------------------------------------------------------------
case(2)
t_c = val1
t_k = t_c + 273.15d0
x_qual = val2
if (t_k < 273.16d0 .or. t_k >= 647.13d0) then
write(*,*) 'ERROR: Saturated state requires T between 0.01 C and 373.9 C'
stop
end if
if (x_qual < 0.0d0 .or. x_qual > 1.0d0) then
write(*,*) 'ERROR: Quality x must be between 0 and 1'
stop
end if
call psat(t_k, p_sat_mpa, rhol, rhov)
p_mpa = p_sat_mpa
p_kpa = p_mpa * 1000.0d0
is_mix = .true.
if (x_qual == 0.0d0) then
state_name = 'Saturated Liquid'
else if (x_qual == 1.0d0) then
state_name = 'Saturated Vapor'
else
state_name = 'Saturated Mixture'
end if
! Calculate properties of saturated liquid
call therm(t_k, rhol, a_term, cjth, cjtt, cp_liq, cv_liq, dpdr, dpdt, &
g_term, h_liq, p_term, s_liq, u_liq)
call viscosity(t_k, rhol, visc_liq)
call thercon(t_k, rhol, therm_liq)
call sound(t_k, p_mpa, sound_liq)
! Calculate properties of saturated vapor
call therm(t_k, rhov, a_term, cjth, cjtt, cp_vap, cv_vap, dpdr, dpdt, &
g_term, h_vap, p_term, s_vap, u_vap)
call viscosity(t_k, rhov, visc_vap)
call thercon(t_k, rhov, therm_vap)
call sound(t_k, p_mpa, sound_vap)
! Blend properties
spec_h = (1.0d0 - x_qual) * h_liq + x_qual * h_vap
spec_s = (1.0d0 - x_qual) * s_liq + x_qual * s_vap
spec_u = (1.0d0 - x_qual) * u_liq + x_qual * u_vap
spec_vol = (1.0d0 - x_qual) / (rhol * 1000.0d0) + x_qual / (rhov * 1000.0d0)
density = 1.0d0 / (spec_vol * 1000.0d0)
cp_val = (1.0d0 - x_qual) * cp_liq + x_qual * cp_vap
cv_val = (1.0d0 - x_qual) * cv_liq + x_qual * cv_vap
speed_sound = (1.0d0 - x_qual) * sound_liq + x_qual * sound_vap
visc_eta = (1.0d0 - x_qual) * visc_liq + x_qual * visc_vap
therm_lambda = (1.0d0 - x_qual) * therm_liq + x_qual * therm_vap
rho = density
! -------------------------------------------------------------
! 3. Pressure and Quality (val1 = P [kPa], val2 = x)
! -------------------------------------------------------------
case(3)
p_kpa = val1
p_mpa = p_kpa / 1000.0d0
x_qual = val2
if (p_mpa < 0.000611d0 .or. p_mpa >= 22.055d0) then
write(*,*) 'ERROR: Saturated state requires P between 0.61 kPa and 22055 kPa'
stop
end if
if (x_qual < 0.0d0 .or. x_qual > 1.0d0) then
write(*,*) 'ERROR: Quality x must be between 0 and 1'
stop
end if
call tsat(p_mpa, t_sat_k, rhol, rhov)
t_k = t_sat_k
t_c = t_k - 273.15d0
is_mix = .true.
if (x_qual == 0.0d0) then
state_name = 'Saturated Liquid'
else if (x_qual == 1.0d0) then
state_name = 'Saturated Vapor'
else
state_name = 'Saturated Mixture'
end if
! Calculate properties of saturated liquid
call therm(t_k, rhol, a_term, cjth, cjtt, cp_liq, cv_liq, dpdr, dpdt, &
g_term, h_liq, p_term, s_liq, u_liq)
call viscosity(t_k, rhol, visc_liq)
call thercon(t_k, rhol, therm_liq)
call sound(t_k, p_mpa, sound_liq)
! Calculate properties of saturated vapor
call therm(t_k, rhov, a_term, cjth, cjtt, cp_vap, cv_vap, dpdr, dpdt, &
g_term, h_vap, p_term, s_vap, u_vap)
call viscosity(t_k, rhov, visc_vap)
call thercon(t_k, rhov, therm_vap)
call sound(t_k, p_mpa, sound_vap)
! Blend properties
spec_h = (1.0d0 - x_qual) * h_liq + x_qual * h_vap
spec_s = (1.0d0 - x_qual) * s_liq + x_qual * s_vap
spec_u = (1.0d0 - x_qual) * u_liq + x_qual * u_vap
spec_vol = (1.0d0 - x_qual) / (rhol * 1000.0d0) + x_qual / (rhov * 1000.0d0)
density = 1.0d0 / (spec_vol * 1000.0d0)
cp_val = (1.0d0 - x_qual) * cp_liq + x_qual * cp_vap
cv_val = (1.0d0 - x_qual) * cv_liq + x_qual * cv_vap
speed_sound = (1.0d0 - x_qual) * sound_liq + x_qual * sound_vap
visc_eta = (1.0d0 - x_qual) * visc_liq + x_qual * visc_vap
therm_lambda = (1.0d0 - x_qual) * therm_liq + x_qual * therm_vap
rho = density
! -------------------------------------------------------------
! 4. Pressure and Enthalpy (val1 = P [kPa], val2 = h [kJ/kg])
! -------------------------------------------------------------
case(4)
p_kpa = val1
p_mpa = p_kpa / 1000.0d0
spec_h = val2
if (p_mpa < 0.000611d0 .or. p_mpa > 100.0d0) then
write(*,*) 'ERROR: Pressure must be between 0.61 kPa and 100000 kPa'
stop
end if
if (p_mpa < 22.055d0) then
! Check if it falls in saturation boundary
call tsat(p_mpa, t_sat_k, rhol, rhov)
call therm(t_sat_k, rhol, a_term, cjth, cjtt, cp_liq, cv_liq, dpdr, dpdt, &
g_term, h_liq, p_term, s_liq, u_liq)
call therm(t_sat_k, rhov, a_term, cjth, cjtt, cp_vap, cv_vap, dpdr, dpdt, &
g_term, h_vap, p_term, s_vap, u_vap)
if (spec_h >= h_liq .and. spec_h <= h_vap) then
! Saturated mixture!
x_qual = (spec_h - h_liq) / (h_vap - h_liq)
t_k = t_sat_k
t_c = t_k - 273.15d0
is_mix = .true.
state_name = 'Saturated Mixture'
call viscosity(t_k, rhol, visc_liq)
call thercon(t_k, rhol, therm_liq)
call sound(t_k, p_mpa, sound_liq)
call viscosity(t_k, rhov, visc_vap)
call thercon(t_k, rhov, therm_vap)
call sound(t_k, p_mpa, sound_vap)
spec_s = (1.0d0 - x_qual) * s_liq + x_qual * s_vap
spec_u = (1.0d0 - x_qual) * u_liq + x_qual * u_vap
spec_vol = (1.0d0 - x_qual) / (rhol * 1000.0d0) + x_qual / (rhov * 1000.0d0)
density = 1.0d0 / (spec_vol * 1000.0d0)
cp_val = (1.0d0 - x_qual) * cp_liq + x_qual * cp_vap
cv_val = (1.0d0 - x_qual) * cv_liq + x_qual * cv_vap
speed_sound = (1.0d0 - x_qual) * sound_liq + x_qual * sound_vap
visc_eta = (1.0d0 - x_qual) * visc_liq + x_qual * visc_vap
therm_lambda = (1.0d0 - x_qual) * therm_liq + x_qual * therm_vap
rho = density
else if (spec_h < h_liq) then
! Subcooled liquid
state_name = 'Subcooled Liquid'
t_min = 273.16d0
t_max = t_sat_k
rho_start = 1.0d0
! Bisection solver for T
do iter = 1, 40
t_mid = 0.5d0 * (t_min + t_max)
call dense(p_mpa, t_mid, rho_start, rho, dpdr)
call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, h_mid, p_term, s_mid, u_term)
if (h_mid < spec_h) then
t_min = t_mid
else
t_max = t_mid
end if
end do
t_k = 0.5d0 * (t_min + t_max)
t_c = t_k - 273.15d0
spec_h = h_mid
spec_s = s_mid
spec_u = u_term
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
else
! Superheated vapor
state_name = 'Superheated Vapor'
t_min = t_sat_k
t_max = 1273.15d0
! Bisection solver for T
do iter = 1, 40
t_mid = 0.5d0 * (t_min + t_max)
rho_start = p_mpa / (gascon() * t_mid)
call dense(p_mpa, t_mid, rho_start, rho, dpdr)
call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, h_mid, p_term, s_mid, u_term)
if (h_mid < spec_h) then
t_min = t_mid
else
t_max = t_mid
end if
end do
t_k = 0.5d0 * (t_min + t_max)
t_c = t_k - 273.15d0
spec_h = h_mid
spec_s = s_mid
spec_u = u_term
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
end if
else
! Supercritical region
state_name = 'Supercritical Fluid'
t_min = 273.16d0
t_max = 1273.15d0
do iter = 1, 40
t_mid = 0.5d0 * (t_min + t_max)
if (t_mid < 647.13d0) then
rho_start = 1.0d0
else
rho_start = p_mpa / (gascon() * t_mid)
end if
call dense(p_mpa, t_mid, rho_start, rho, dpdr)
call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, h_mid, p_term, s_mid, u_term)
if (h_mid < spec_h) then
t_min = t_mid
else
t_max = t_mid
end if
end do
t_k = 0.5d0 * (t_min + t_max)
t_c = t_k - 273.15d0
spec_h = h_mid
spec_s = s_mid
spec_u = u_term
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
end if
! -------------------------------------------------------------
! 5. Pressure and Entropy (val1 = P [kPa], val2 = s [kJ/(kg K)])
! -------------------------------------------------------------
case(5)
p_kpa = val1
p_mpa = p_kpa / 1000.0d0
spec_s = val2
if (p_mpa < 0.000611d0 .or. p_mpa > 100.0d0) then
write(*,*) 'ERROR: Pressure must be between 0.61 kPa and 100000 kPa'
stop
end if
if (p_mpa < 22.055d0) then
! Check if it falls in saturation boundary
call tsat(p_mpa, t_sat_k, rhol, rhov)
call therm(t_sat_k, rhol, a_term, cjth, cjtt, cp_liq, cv_liq, dpdr, dpdt, &
g_term, h_liq, p_term, s_liq, u_liq)
call therm(t_sat_k, rhov, a_term, cjth, cjtt, cp_vap, cv_vap, dpdr, dpdt, &
g_term, h_vap, p_term, s_vap, u_vap)
if (spec_s >= s_liq .and. spec_s <= s_vap) then
! Saturated mixture!
x_qual = (spec_s - s_liq) / (s_vap - s_liq)
t_k = t_sat_k
t_c = t_k - 273.15d0
is_mix = .true.
state_name = 'Saturated Mixture'
call viscosity(t_k, rhol, visc_liq)
call thercon(t_k, rhol, therm_liq)
call sound(t_k, p_mpa, sound_liq)
call viscosity(t_k, rhov, visc_vap)
call thercon(t_k, rhov, therm_vap)
call sound(t_k, p_mpa, sound_vap)
spec_h = (1.0d0 - x_qual) * h_liq + x_qual * h_vap
spec_u = (1.0d0 - x_qual) * u_liq + x_qual * u_vap
spec_vol = (1.0d0 - x_qual) / (rhol * 1000.0d0) + x_qual / (rhov * 1000.0d0)
density = 1.0d0 / (spec_vol * 1000.0d0)
cp_val = (1.0d0 - x_qual) * cp_liq + x_qual * cp_vap
cv_val = (1.0d0 - x_qual) * cv_liq + x_qual * cv_vap
speed_sound = (1.0d0 - x_qual) * sound_liq + x_qual * sound_vap
visc_eta = (1.0d0 - x_qual) * visc_liq + x_qual * visc_vap
therm_lambda = (1.0d0 - x_qual) * therm_liq + x_qual * therm_vap
rho = density
else if (spec_s < s_liq) then
! Subcooled liquid
state_name = 'Subcooled Liquid'
t_min = 273.16d0
t_max = t_sat_k
rho_start = 1.0d0
! Bisection solver for T
do iter = 1, 40
t_mid = 0.5d0 * (t_min + t_max)
call dense(p_mpa, t_mid, rho_start, rho, dpdr)
call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, h_mid, p_term, s_mid, u_term)
if (s_mid < spec_s) then
t_min = t_mid
else
t_max = t_mid
end if
end do
t_k = 0.5d0 * (t_min + t_max)
t_c = t_k - 273.15d0
spec_h = h_mid
spec_s = s_mid
spec_u = u_term
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
else
! Superheated vapor
state_name = 'Superheated Vapor'
t_min = t_sat_k
t_max = 1273.15d0
! Bisection solver for T
do iter = 1, 40
t_mid = 0.5d0 * (t_min + t_max)
rho_start = p_mpa / (gascon() * t_mid)
call dense(p_mpa, t_mid, rho_start, rho, dpdr)
call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, h_mid, p_term, s_mid, u_term)
if (s_mid < spec_s) then
t_min = t_mid
else
t_max = t_mid
end if
end do
t_k = 0.5d0 * (t_min + t_max)
t_c = t_k - 273.15d0
spec_h = h_mid
spec_s = s_mid
spec_u = u_term
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
end if
else
! Supercritical region
state_name = 'Supercritical Fluid'
t_min = 273.16d0
t_max = 1273.15d0
do iter = 1, 40
t_mid = 0.5d0 * (t_min + t_max)
if (t_mid < 647.13d0) then
rho_start = 1.0d0
else
rho_start = p_mpa / (gascon() * t_mid)
end if
call dense(p_mpa, t_mid, rho_start, rho, dpdr)
call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, h_mid, p_term, s_mid, u_term)
if (s_mid < spec_s) then
t_min = t_mid
else
t_max = t_mid
end if
end do
t_k = 0.5d0 * (t_min + t_max)
t_c = t_k - 273.15d0
spec_h = h_mid
spec_s = s_mid
spec_u = u_term
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
end if
end select
! Density conversion: rho is in g/cm3. Output specific volume v in m3/kg.
if (.not. is_mix) then
density = rho * 1000.0d0 ! kg/m3
spec_vol = 1.0d0 / density ! m3/kg
end if
! If it's single phase, compute dynamic viscosity in Pa s
! visc_eta is in microPascal-seconds. Let's convert to Pa s: * 1.0E-6
visc_eta = visc_eta * 1.0d-6
! Convert thermal conductivity to W/(m K): therm_lambda is in mW/(m K) => / 1000.d0
therm_lambda = therm_lambda / 1000.d0
! Convert quality
if (.not. is_mix) then
if (state_name == 'Subcooled Liquid') then
x_qual = 0.0d0
else if (state_name == 'Superheated Vapor') then
x_qual = 1.0d0
else
x_qual = -1.0d0 ! Supercritical has no quality defined
end if
end if
! Print output in parsable format
write(*,'(A,A)') 'STATE=', trim(state_name)
write(*,'(A,F14.4)') 'T=', t_c
write(*,'(A,F14.4)') 'P=', p_kpa
write(*,'(A,F14.6)') 'V=', spec_vol
write(*,'(A,F14.4)') 'H=', spec_h
write(*,'(A,F14.6)') 'S=', spec_s
write(*,'(A,F14.4)') 'U=', spec_u
write(*,'(A,F14.4)') 'CP=', cp_val
write(*,'(A,F14.4)') 'CV=', cv_val
write(*,'(A,F14.2)') 'SOUND=', speed_sound
write(*,'(A,E14.6)') 'VISC=', visc_eta
write(*,'(A,F14.6)') 'THERM=', therm_lambda
write(*,'(A,F14.4)') 'QUALITY=', x_qual
end program steam_calc
Solver Description
This high-precision solver implements the 1984 NBS/NRC formulations (Haar, Gallagher, Kell) for the thermodynamic properties of water. It evaluates the Helmholtz free energy function $A(\rho, T)$ and its analytical derivatives to calculate density, specific volume, enthalpy, entropy, internal energy, speed of sound, specific heats ($C_p, C_v$), viscosity, and thermal conductivity. An iterative numerical bisection solver is used in the backend to solve for temperature ($T$) and density ($\rho$) when enthalpy ($h$) or entropy ($s$) are provided as input parameters.
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
🛠️ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -ffree-line-length-none steam.f90 steam_calc.f90 -o steam_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
steam_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Input combination mode (1=T/P, 2=T/x, 3=P/x, 4=P/h, 5=P/s)
1
! Parameter 1 (Temperature [°C] or Pressure [kPa])
150.0
! Parameter 2 (Pressure [kPa] or Quality [0-1] or Enthalpy [kJ/kg] or Entropy [kJ/kg-K])
101.325
1
! Parameter 1 (Temperature [°C] or Pressure [kPa])
150.0
! Parameter 2 (Pressure [kPa] or Quality [0-1] or Enthalpy [kJ/kg] or Entropy [kJ/kg-K])
101.325