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