diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..7616378 --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +build/ +fortran/*.for +mixtures/*.MIX +fluids/*.FLD +fluids/*.PPF +fluids/HMX.BNC +fortran/MANUAL.TXT diff --git a/CMakeLists.txt b/CMakeLists.txt index cf0e901..64a39a4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,65 +1,98 @@ -cmake_minimum_required(VERSION 2.8) -project(REFPROP Fortran) - -SET (USE_OPENMP ON CACHE BOOL "Use OpenMP") -SET (CMAKE_BUILD_TYPE "RELEASE" CACHE INTERNAL "Set test build") - -SET(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake/Modules/") -IF(NOT CMAKE_Fortran_COMPILER_SUPPORTS_F90) - MESSAGE(FATAL_ERROR "Fortran compiler does not support F90") -ENDIF(NOT CMAKE_Fortran_COMPILER_SUPPORTS_F90) +cmake_minimum_required(VERSION 3.3) +SET(CMAKE_BUILD_TYPE Debug) -INCLUDE(${CMAKE_MODULE_PATH}/SetFortranFlags.cmake) -INCLUDE(${CMAKE_MODULE_PATH}/SetParallelizationLibrary.cmake) -find_package(PythonLibs 3 REQUIRED) -find_package(PythonInterp) +project(RGPTableGen Fortran) -# running on windows and refprop is installed? then copy what we need -if (WIN32) - # Set locations to store output - set(REFPROP_FLUID_FOLDER "${CMAKE_BINARY_DIR}") - set(REFPROP_DLL_FILE "${CMAKE_BINARY_DIR}") - - # Refprop (new version) is a 64 bit binary - find_path(REFPROP_BASE_PATH "REFPRP64.dll" "c:/program files/REFPROP" "c:/program files (x86)/REFPROP") +if (NOT CMAKE_ARCHIVE_OUTPUT_DIRECTORY) + set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) +endif() + +if (NOT CMAKE_LIBRARY_OUTPUT_DIRECTORY) + set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) +endif() + +if (NOT CMAKE_RUNTIME_OUTPUT_DIRECTORY) + set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin) endif() -if (USE_OPENMP) - message(STATUS "REFPROP: Each REFPROP instance may use up to ${OMP_NUM_PROCS} threads with OpenMP") +# Build Environment +SET(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules/") + +if (NOT REFPROP_FLAGS_ADDED) + set(REFPROP_FLAGS_ADDED "true" CACHE INTERNAL "Has refprop flags been added?") + + if(NOT CMAKE_Fortran_COMPILER_SUPPORTS_F90) + MESSAGE(FATAL_ERROR "Fortran compiler does not support F90") + else() + SET (CMAKE_BUILD_TYPE "DEBUG" CACHE INTERNAL "Set test build") + INCLUDE("${CMAKE_MODULE_PATH}/SetFortranFlags.cmake") + INCLUDE("${CMAKE_MODULE_PATH}/SetCompileFlag.cmake") + endif() + + SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}" Fortran "-fpic") + SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}" Fortran "-ffast-math") + SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}" Fortran "-fno-common") + SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}" Fortran "-ffloat-store") + SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}" Fortran "-fall-intrinsics") + + message(STATUS "${CMAKE_Fortran_FLAGS_RELEASE}") endif() -# Copy fluid folder to output directory -if (REFPROP_BASE_PATH) - message(STATUS "REFPROP: REFPROP found at ${REFPROP_BASE_PATH}, copying/using existing fluid/fortran files") -else() - set(REFPROP_BASE_PATH "${CMAKE_CURRENT_SOURCE_DIR}" CACHE INTERNAL "Code Path") -endif() -# Remove non required files -message("${REFPROP_BASE_PATH}") +# Options +SET (USE_OPENMP OFF CACHE BOOL "Use OpenMP") + +find_package(PythonInterp REQUIRED) + +##################### +# CONFIGURE REFPROP # +##################### + +set(REFPROP_BASE_PATH "${CMAKE_CURRENT_SOURCE_DIR}" CACHE INTERNAL "Refprop code location") + +# Set locations to store output +set(REFPROP_FLUID_FOLDER "${CMAKE_BINARY_DIR}") +set(REFPROP_DLL_FILE "${CMAKE_BINARY_DIR}") + + +# Copy REFPROP fluid files +message(STATUS "REFPROP: REFPROP fluids found at ${REFPROP_BASE_PATH}, copying fluid files") file(COPY "${REFPROP_BASE_PATH}/fluids" DESTINATION "${CMAKE_BINARY_DIR}") +message(STATUS "REFPROP: Fluid files copied") + +# Get fortran sources FILE(GLOB refprop_files "${REFPROP_BASE_PATH}/fortran/*.FOR") + +# Remove non-required items LIST(REMOVE_ITEM refprop_files -"${REFPROP_BASE_PATH}/fortran/COMMONS.FOR" -"${REFPROP_BASE_PATH}/fortran/COMTRN.FOR" +"${REFPROP_BASE_PATH}/fortran/COMMONS.FOR" # include file +"${REFPROP_BASE_PATH}/fortran/COMTRN.FOR" # include file "${REFPROP_BASE_PATH}/fortran/UTILITY.FOR" ) -execute_process(COMMAND ${PYTHON_EXECUTABLE} "patchUtility.py" "${REFPROP_BASE_PATH}/fortran/UTILITY.FOR" "${CMAKE_CURRENT_SOURCE_DIR}/UTILITY.FOR") -LIST(APPEND refprop_files "${CMAKE_CURRENT_SOURCE_DIR}/UTILITY.FOR") -message("${refprop_files}") +# Make a copy of files (which are to be patched) +message(STATUS "REFPROP: Patching fortran files") +file(COPY "${REFPROP_BASE_PATH}/fortran/COMMONS.FOR" DESTINATION "${CMAKE_BINARY_DIR}/fortran") +file(COPY "${REFPROP_BASE_PATH}/fortran/COMTRN.FOR" DESTINATION "${CMAKE_BINARY_DIR}/fortran") +file(RENAME "${CMAKE_BINARY_DIR}/fortran/COMMONS.FOR" "${CMAKE_BINARY_DIR}/fortran/commons.for") +file(RENAME "${CMAKE_BINARY_DIR}/fortran/COMTRN.FOR" "${CMAKE_BINARY_DIR}/fortran/comtrn.for") -# Add some nice flags -if(UNIX OR MINGW) - SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}" Fortran "${CMAKE_Fortran_FLAGS} -ffast-math -fno-common -ffloat-store" FORCE) -endif() +execute_process(COMMAND ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/patchUtility.py" + "${REFPROP_BASE_PATH}/fortran/UTILITY.FOR" + "${CMAKE_BINARY_DIR}/UTILITY.FOR") + +LIST(APPEND refprop_files "${CMAKE_BINARY_DIR}/UTILITY.FOR") + +# Set include directory +include_directories("${CMAKE_BINARY_DIR}/fortran") # Setup build include_directories("${REFPROP_BASE_PATH}/fortran") # Use common core method add_library(REFPROPObjects OBJECT ${refprop_files}) -add_library(REFPRP64 SHARED $) -add_executable(RGPTableGen $ "${CMAKE_CURRENT_SOURCE_DIR}/RGPgenerator.F90") +add_executable(RGPTableGen $ "${CMAKE_CURRENT_SOURCE_DIR}/RGP.for") + +target_compile_options(RGPTableGen PRIVATE -std=f2008 -ffree-form) diff --git a/README.md b/README.md index 48edd0e..138a076 100644 --- a/README.md +++ b/README.md @@ -3,20 +3,58 @@ RGP Generator This program was developed to generate custom step size RGP files using the NIST fluid properties. -Included Files --------------- -fluids(folder) - Includes REFPROP fluids needed for generating the rgp file -Source Code(folder) - Includes the source code for the program (within RGPgenerator.F90, the rest are REFPROP subroutines, and RGP_Generator.prj is the simple fortran project file) -CFX Solver Guide.pdf - ANSYS CFX Documentation which was useful when formatting the RGP gen output file (See section 12.6) -RGP_Generator64.exe - 64 Bit version of the program -RGP_Generator32.exe - 32 Bit version of the program -README.txt - -How to use ------------ -1) Open either the 32 or 64 bit version of RGP_Generator.exe -2) Follow the steps to generate the output file -3) Upload the RGP file to ANSYS CFX +How to Use +---- + +RGP File gen is a simple program generating complex files for use with the CFX fluid dynamics package. +Once the program is built (see below about how to build) you may run the program on the command line +and answer the questions with regard to your current project. Be advised that the temperature and pressure +ranges that you enter must cover both the total and static states expected in the whole flow field else +you will end up with range clipping (as will be seen in the CFX log file). + +New in this version is the ability to put the answers to the questions in the program arguments. Please be +aware that the order matters, and the order is the same as the normal question asking order. You will see +the output however printed in the terminal so you will be made aware of any errors. + +How to Build +===== + +* Bad news: Refprop is licenced by NIST, and should not be distributed. +* Good news: We have everything sorted for you! + +RGP Generator is simple to build, your computer may even have everything required +already installed. Please use the following steps to get setup. Be aware if you have +linux, that you would have to had REFPROP installed on windows to get the code and fluid +files. + +1) Copy all fluids from your REFPROP install path (c:/program files/refprop) to the fluids folder +2) Copy all the fortran files to the fortran directory +3) Follow the steps below based on your OS + +How to Build (Windows) +----- +1) Install mingw with fortran, if this is forign then watch the video available [here](https://www.youtube.com/watch?v=oVfAU1ziOjg +2) Install cmake from [here](https://cmake.org/download/) +3) Make a folder called build in the current directory +4) Run cmake GUI +* Set `where is the source code located` to the folder with RGP.for +* Set `where to build the binaries` to the build folder you made +* Press configure and select `MingGW Makefiles` in `Specify the generator for this project` +* Press generate +5) Open a new command window and ensure the current directory is the build folder +6) Run `mingw32-make` which should be located in the bin directory of where you installed mingw32 + +How to Build (Linux) +------- + +Linux is much nicer! + +1) Install fortran (`apt-get install gfortran` or `pacman -S gfortran`) based on your system +2) Install cmake (again using apt or pacman) +3) Navigate to the source folder and make a build directory +4) In the build directory run cmake .. +5) Run make + Licence ======= diff --git a/RGP.for b/RGP.for new file mode 100644 index 0000000..2be4621 --- /dev/null +++ b/RGP.for @@ -0,0 +1,903 @@ +! This program was written to generate Real Gas Property (RGP) files to be used +! with ANSYS CFX (Current version 14.5). The user is required to set certain +! input variables which generates an RGP file that follows the ANSYS guidelines +! +! See Anysys manual section "12.7.1.1. Detailed .rgp File Format" +! this is easily found on Google. + +! Current Features: +! Output file to ANSYS CFX version 14.5 specifications +! Input using arguments of interactively +! Any fluid :D + +! Notes: +! - Must be compiled with -std=f2008 -ffree-form +! - Other flags can be whatever +! - Must be linked to refprop objects + +module RGP_DATA + integer ncmax + parameter (ncmax=20) !max number of components in mixture + double precision, dimension (ncmax) :: x(ncmax),xliq(ncmax),xvap(ncmax) + double precision, allocatable :: t(:),p(:),h(:,:),w(:,:),v(:,:),cv(:,:),cp(:,:),s(:,:),dPdv(:,:),eta(:,:),tcx(:,:), & + &tsat(:),hsat(:),wsat(:),vsat(:),cvsat(:),cpsat(:),ssat(:),dPdvsat(:),etasat(:),tcxsat(:), & + &tsat2(:),hsatl(:),wsatl(:),rhosatl(:),cvsatl(:),cpsatl(:),ssatl(:),dPdrhosatl(:),etasatl(:),tcxsatl(:), & + &psat2(:),hsatv(:),wsatv(:),rhosatv(:),cvsatv(:),cpsatv(:),ssatv(:),dPdrhosatv(:),etasatv(:),tcxsatv(:), & + &blank1(:),blank2(:),arr(:) + double precision :: acf,conp,conv,conpl,convl,conpv,convv,cpcrit,cvcrit,density,dpdrho,q,rgas,rhocrit,rhol,rhov, & + &dc,dip,dl,dm,dmax,dpdrhol,dpdrhov,dv,e,enthalpy,enthalpyl,enthalpyv,pstep,pstepsat,scrit,speedcrit, & + &entropy,entropyl,entropyv,hcrit,htyp,pc,pm,pmax,pmaxsat,pminsat,psat,preslower,upperPres,speedOsound, & + &speedOsoundl,speedOsoundv,supercool,tc,templower,uppertemp,TminSat,TmaxSat,thermalconductivity,thermalconductivityl, & + &thermalconductivityv,thermcrit,tliq,tm,tmax,tmin,tnbp,tpp,ttp,tstep,tvap,visccrit,viscosity,viscosityl,viscosityv,wm, & + &y,zc,dbl + integer :: i,j,ierr,kq,np,nt,nsat,npsat,ilng + character :: hrf*3, herr*255, lfluid*30, loutfile*30 + character :: lFileName*255, hnam*12, hn80*80, hstr*255 + character :: hcasn*12 + character(LEN=255) :: hf(ncmax) + character(LEN=7) :: hfmix + integer :: prevStatus + integer :: currStatus + integer :: tenPercent + integer :: AllocStatus +end module + +program RGP + USE RGP_DATA + implicit none + + write(*,*) ' Carleton University' + write(*,*) ' Department of Mechanical and Aerospace Engineering' + write(*,*) ' &' + write(*,*) ' Queensland University of Technology' + write(*,*) ' Department of Chemistry Physics and Mechanical Engineering' + write(*,*) + write(*,*) ' 88888888ba ,ad8888ba, 88888888ba ' + write(*,*) ' 88 "8b d8" "8b 88 "8b ' + write(*,*) ' 88 .8P d8 88 8P ' + write(*,*) ' 88aaaaaa8P 88 88aaaaaa8P ' + write(*,*) ' 88""""88 88 88888 88"""""" ' + write(*,*) ' 88 `8b Y8, `88 88 ' + write(*,*) ' 88 `8b Y8a. .a88 88 ' + write(*,*) ' 88 `8b `"Y88888P" 88 ' + write(*,*) + write(*,*) ' FILE GENERATOR' + write(*,*) + + call checkArg(); + call getFilesConfig(); + call setupREFPROP(); + call getRGPBounds(); + + write(*,*)'Writing header...' + call writeRGPHeader(); + + write(*,*)'Allocating arrays...' + call allocateArrays(); + + write(*,*)'Calc sat table...' + call calculateSatTable(); + + write(*,*)'Calc sat table...' + call calculateSuperTable(); + + write(*,*)'Writing super table...' + call writeSuperTable(); + + write(*,*)'Writing sat table...' + call writeSatTable(); + close(10); + +contains + subroutine checkArg() + + if (COMMAND_ARGUMENT_COUNT() .gt. 1) then + write(*,*) 'Using arguments for inputs' + call get_command_argument(1, lOutfile) + call get_command_argument(2, lFileName) + end if + + if (COMMAND_ARGUMENT_COUNT() .eq. 2+9) then + call getDblArg(3, TempLower) + call getDblArg(4, UpperTemp) + call getIntArg(5, nt) + call getDblArg(6, PresLower) + call getDblArg(7, UpperPres) + call getIntArg(8, np) + call getDblArg(9, TminSat) + call getDblArg(10, TmaxSat) + call getIntArg(11, npsat) + end if + + end subroutine checkArg + + subroutine getFilesConfig() + USE RGP_DATA + implicit none + + ! Check if lOutfils is set, or get from cmd line + i = SCAN(lOutfile,'.') + if (i == 0) then + write(*,'(A)') 'Output Data File Name (i.e. filename.txt):' + read (*,*) lOutfile + else + write(*,*) 'Output Data File: ', lOutfile + end if + + ! Check if lFileName is set, or get from cmd line + i = SCAN(lFileName,'.') + if (i == 0) then + write(*,'(A)') 'Fluid file name (i.e. filename.fld):' + read (*,*) lFileName + else + write(*,*) 'Fluid file name: ', lFileName + end if + + ! Trim excess spaces + lOutfile = TRIM(lOutfile) + lFileName = TRIM(lFileName) + + ! Set lfluid from the file name + i = SCAN(lFileName,'.') + lfluid = lFileName(1:i-1) + + ! Create write unit for output file + OPEN(UNIT=10,FILE=lOutFile) + end subroutine getFilesConfig + + subroutine setupREFPROP() + USE RGP_DATA + implicit none + + ! Set fluid to load + i=1 + hf(1)=lFileName + + ! default fluid path is fluids/ + call SETPATH ('fluids/') + + ! Use default mixtures setup and + ! fluids preferred (DEFault) EOS + hfmix='hmx.bnc' + hrf='DEF' + + ! Call refprop setup + call SETUP (i,hf,hfmix,hrf,ierr,herr) + call checkRPError(ierr, herr) + + write(*,*) 'REFPROP loaded' + + call INFO (1,wm,ttp,tnbp,tc,pc,dc,zc,acf,dip,rgas) + call NAME (1,hnam,hn80,hcasn) + call PASSCMN ('ptpn',0,1,0,hstr,ilng,tpp,arr,ierr,herr) + call MAXT (x,tm,pm,Dm,ierr,herr) + call LIMITS(htyp,x,tmin,tmax,Dmax,pmax) + + dc = wm * dc + Dmax = wm * Dmax + hnam = TRIM(hnam); + + write(*,*) + write(*,*) 'Short name: ',hnam + write(*,*) 'Molecular Weight: ',wm,'g/mol' + write(*,*) 'Triple Point Temperature: ',ttp,'K' + write(*,*) 'Triple Point Pressure: ',tpp,'kPa' + write(*,*) 'Normal Boiling Point: ',tnbp,'K' + + write(*,*) 'Critical Temperature: ',tc,'K' + write(*,*) 'Critical Pressure: ',pc,'kPa' + write(*,*) 'Critical Density: ',dc,'kg/m^3' + write(*,*) 'Compressibility at Critical Point: ',zc + write(*,*) 'Gas Constant: ',rgas,'J/mol-K' + + write(*,*) + write(*,*)'Limits for this fluid:' + write(*,*)'Min Temperature: ',tmin,'K' + write(*,*)'Max Temperature: ',tmax,'K' + write(*,*)'Max Saturation temperature: ',tm,'K' + write(*,*)'Max Density: ',Dmax,'kg/m^3' + write(*,*)'Max Pressure: ',pmax,'kPa' + end subroutine setupREFPROP + + subroutine getRGPBounds() + USE RGP_DATA + implicit none + + ! PROMPT FOR BOUNDS AND STEP SIZE + write(*,*)'--------------------------------------------' + write(*,*) + call askQuestion('Lower bound for Temperature [K]:', TempLower) + call askQuestion('Upper bound for Temperature [K]:', UpperTemp) + call askQuestion_int('Number of Steps for Temperature:', nt) + + TempLower = max(TempLower, tmin); + UpperTemp = min(UpperTemp, tmax); + + call askQuestion('Lower bound for Pressure [Pa]:', PresLower) + call askQuestion('Upper bound for Pressure [Pa]:', UpperPres) + call askQuestion_int('Number of Steps for Pressure :', np) + + PresLower = max(PresLower, DBLE(0.0)); + UpperPres = min(UpperPres, pmax*1000); + + call askQuestion('Minimum Temperature for Saturation Tables [k]:', TminSat) + call askQuestion('Maximum Temperature for Saturation Tables [k]:', TmaxSat) + call askQuestion_int('Number of Steps for Saturation Tables :', npsat) + + TminSat = max(TminSat, tmin); + TmaxSat = min(TmaxSat, tm); + + ! set t and p step + tstep = (UpperTemp-TempLower)/REAL(nt-1) + pstep = (UpperPres-PresLower)/REAL(np-1) + end subroutine getRGPBounds + + !---------------------------------------------------------------------! + ! DATASET CALCULATION ! + ! ! + ! getSatLineData : Gets information on the saturation line ! + ! getSuperDataTP : Gets information at any vapour point ! + ! getSuperDataPQ : calls getSatLineData for super table ! + ! ! + !---------------------------------------------------------------------! + + + subroutine calculateSuperTable() + USE RGP_DATA + implicit none + + prevStatus = 0 + currStatus = 0 + tenPercent = int(UpperPres-PresLower)/10 + j=1 + p(j) = PresLower + do while (j <= np) + i=1 + t(i) = TempLower + + ! show progress every 10% + currStatus = int(p(j) - PresLower)/tenPercent + if (prevStatus /= currStatus) then + write (*,*) (currStatus*10), '%' + prevStatus = currStatus + end if + + do while (i <= nt) + + ! if p(j) <= pc (below critical pressure), calc Phi saturated line + if (p(j) <= pc*1000) then + PSat = p(j)/1000 + q = 1 + + ! find the sat temp for this pressure + call PQFLSH (PSat,q,x,kq,tvap,rhov,Dl,Dv,xliq,xvap,e,enthalpyv,entropyv,convv,conpv,speedOsoundv,ierr,herr) + call checkRPError(ierr, herr) + + ! if t(i) <= tvap (below sat temp), use phi(Tsat, p), q=1 line + if (t(i) <= tvap) then + call getSuperDataPQ(PSat, h(i,j), w(i,j), v(i,j), cv(i,j), cp(i,j),& + & s(i,j), dPdv(i,j),eta(i,j), tcx(i,j)) + + ! if t(i) > tvap (above sat temp), get vapour properties + else + call getSuperDataTP(t(i), p(j), h(i,j), w(i,j), v(i,j), cv(i,j), cp(i,j),& + & s(i,j), dPdv(i,j), eta(i,j), tcx(i,j)) + end if + + ! if p(j) > pc (above critical pressure), get vapour properties + else + call getSuperDataTP(t(i), p(j), h(i,j), w(i,j), v(i,j), cv(i,j), cp(i,j),& + & s(i,j), dPdv(i,j),eta(i,j), tcx(i,j)) + end if + + i = i + 1 + if (i <= nt) then + t(i) = t(i-1)+tstep + end if + end do + + ! Fill 1D Table (Sat line) + if (p(j)<(Pc*1000))then + q = 1 + PSat = p(j)/1000 + call getSuperDataPQ(PSat, hsat(j), wsat(j), vsat(j), cvsat(j), cpsat(j),& + & ssat(j), dPdvsat(j),etasat(j), tcxsat(j), tsat(j)) + else + ! fill with phi(Tmin, P) + q=1 + PSat = p(j)/1000 + tsat(j) = TempLower + call getSuperDataTP(TempLower, p(j), hsat(j), wsat(j), vsat(j), cvsat(j), cpsat(j),& + & ssat(j), dPdvsat(j),etasat(j), tcxsat(j)) + end if + + j = j+1 + if (j <= np) then + p(j) = p(j-1)+pstep + end if + end do + + end subroutine calculateSuperTable + + + subroutine calculateSatTable() + write (*,*) 'Filling saturation table...' + + ! Get pressure step + q=1 + call TQFLSH (TminSat,q,x,kq,PminSat,density,Dl,Dv,x,y,e,enthalpy,entropy,conv,conp,speedOsound,ierr,herr) + call checkRPError(ierr, herr) + + call TQFLSH (TmaxSat,q,x,kq,PmaxSat,density,Dl,Dv,x,y,e,enthalpy,entropy,conv,conp,speedOsound,ierr,herr) + call checkRPError(ierr, herr) + + pstepsat = (Pmaxsat-Pminsat)/REAL(npsat-1) + + ! Iterate + j=1 + psat2(j) = PminSat + do while (j <= npsat) + + q=0 + call getSatLineData(psat2(j), hsatl(j), wsatl(j), rhosatl(j), cvsatl(j), cpsatl(j),& + & ssatl(j), dPdrhosatl(j), etasatl(j), tcxsatl(j)) + + q=1 + call getSatLineData(psat2(j), hsatv(j), wsatv(j), rhosatv(j), cvsatv(j), cpsatv(j),& + & ssatv(j), dPdrhosatv(j), etasatv(j), tcxsatv(j), tsat2(j)) + + blank1(j)=0 + blank2(j)=0 + + j = j+1 + if (j <= npsat) then + psat2(j) = psat2(j-1) + pstepsat + end if + end do + end subroutine calculateSatTable + + !---------------------------------------------------------------------! + ! FLID STATE CALCULATION ! + ! ! + ! getSatLineData : Gets information on the saturation line ! + ! getSuperDataTP : Gets information at any vapour point ! + ! getSuperDataPQ : calls getSatLineData for super table ! + ! ! + !---------------------------------------------------------------------! + + subroutine getSatLineData(psat2_j, hsat_q, wsat_q, rhosat_q, cvsat_q, cpsat_q, ssat_q, dPdrhosat_q, etasat_q, tcxsat_q, tsat_q) + USE RGP_DATA + implicit none + + double precision, intent(in) :: psat2_j + double precision, intent(out) :: hsat_q, wsat_q, rhosat_q, cvsat_q, cpsat_q, ssat_q, dPdrhosat_q, etasat_q, tcxsat_q + double precision, optional, intent(out) :: tsat_q + double precision :: rho + + call PQFLSH (psat2_j,q,x,kq,tvap,rho,Dl,Dv,xliq,xvap,e,enthalpy,entropy,conv,conp,speedOsound,ierr,herr) + call checkRPError(ierr, herr) + + call DPDD (tvap,rho,x,dPdrho) + call checkRPError(ierr, herr) + + call TRNPRP (tvap,rho,x,viscosity,thermalconductivity,ierr,herr) + call checkRPError(ierr, herr) + + rhosat_q = rho*wm + hsat_q = enthalpy*1000/wm + ssat_q = entropy*1000/wm + cvsat_q = conv*1000/wm + cpsat_q = conp*1000/wm + wsat_q = speedOsound + + dPdrhosat_q = wm/(dPdrho*1000) + etasat_q = viscosity/1000000 + tcxsat_q = thermalconductivity + + if (present(tsat_q)) tsat_q = tvap + end subroutine getSatLineData + + + subroutine getSuperDataPQ(p_s, h_s, w_s, v_s, cv_s, cp_s, s_s, dPdv_s, eta_s, tcx_s, tvap_s) + USE RGP_DATA + implicit none + + double precision, intent(in) :: p_s + double precision, intent(out) :: h_s, w_s, v_s, cv_s, cp_s, s_s, dPdv_s, eta_s, tcx_s + double precision, optional, intent(out) :: tvap_s + double precision :: rho, dPdrhosat_q + + call getSatLineData(p_s, h_s, w_s, rho, cv_s, cp_s, s_s, dPdrhosat_q, eta_s, tcx_s, tvap_s) + + v_s = 1/rho + dPdv_s = -((rho*wm)**2)/dPdrhosat_q ! (dPdrho/wm*1000*(-(rho*wm)**2)) + end subroutine getSuperDataPQ + + + subroutine getSuperDataTP(t_s, p_s, h_s, w_s, v_s, cv_s, cp_s, s_s, dPdv_s, eta_s, tcx_s) + USE RGP_DATA + implicit none + + double precision, intent(in) :: t_s, p_s + double precision, intent(out) :: h_s, w_s, v_s, cv_s, cp_s, s_s, dPdv_s, eta_s, tcx_s + double precision :: rho + + call TPFLSH (t_s,p_s/1000,x,rho,dl,dv,xliq,xvap,q,e,enthalpy,entropy,conv,conp,speedOsound,ierr,herr) + call checkRPError(ierr, herr) + + call DPDD (t_s,rho,x,dPdrho) + call checkRPError(ierr, herr) + + call TRNPRP (t_s,rho,x,viscosity,thermalconductivity,ierr,herr) + call checkRPError(ierr, herr) + + h_s = enthalpy*1000/wm + w_s = speedOsound + v_s = 1/(rho*wm) + cv_s = conv*1000/wm + cp_s = conp*1000/wm + s_s = entropy*1000/wm + dPdv_s = (dPdrho/wm*1000*(-(rho*wm)**2)) + eta_s = viscosity/1000000 + tcx_s = thermalconductivity + end subroutine getSuperDataTP + + + !---------------------------------------------------------------------! + ! TABLE WRITING FUNCTIONS ! + ! ! + ! writeRGPHeader : Writes RGP main and sub header ! + ! ! + ! See below for other table writing functions ! + !---------------------------------------------------------------------! + + + subroutine writeRGPHeader() + USE RGP_DATA + implicit none + + ! SETUP outPUT FILE + write(10,901)'$$$$HEADER' + + do i=1,2 + if (i .eq. 2) write(10,901)'$$$$DATA' + + write(10,901)'$$$'//TRIM(lfluid) !character*8 (key into the RGP file see NOTE 3) if this program is used with a fluid other than CO2, this line will need to be changed, along with the corresponding line in the $$$DATA header + write(10,*) 1 !enter an integer (this line is ignored in ANSYS CFX) + write(10,901)'$$PARAM' ! See NOTE 4 + write(10,*) 26 !enter an integer(number of parameters) + write(10,901)'DESCRIPTION' + write(10,901)TRIM(lfluid)//' from NIST'!character*50(description of the material) + write(10,901)'NAME' + write(10,901)TRIM(lfluid) !character*8(material name, same as $$$) + write(10,901)'inDEX' + write(10,901)TRIM(lfluid) !character*50 (index into clients RGDB program) + write(10,901)'MODEL' + write(10,*)3 !integer (level of property info available 1,2,or 3) + write(10,901)'UNITS' + write(10,*)1 !integer(unit system of 1 ,2,3,4 or 5) + write(10,901)'PMin_SUPERHEAT' ! at Triple Point Pressure? + write(10,900)PresLower !real (Pa) + write(10,901)'PMAX_SUPERHEAT' ! at Critical Pressure? + write(10,900)UpperPres !real (Pa) + write(10,901)'TMin_SUPERHEAT' ! at Triple Point Temperature? + write(10,900)TempLower !real (K) + write(10,901)'TMAX_SUPERHEAT' ! at Critical Temperature? + write(10,900)UpperTemp !real (K) + write(10,901)'TMin_SATURATION' + write(10,900)TminSat + write(10,901)'TMAX_SATURATION' + write(10,900)TmaxSat + write(10,901)'SUPERCOOLING' + write(10,900)SuperCool !real(supercooling level in superheat tables)(Optional) + write(10,901)'P_CRITICAL' + write(10,900)Pc*1000 !real + write(10,901)'P_TRIPLE' ! See NOTE 7 + write(10,900)tpp*1000!real + write(10,901)'T_CRITICAL' + write(10,900)Tc!real + write(10,901)'T_TRIPLE' + write(10,900)Ttp!real + write(10,901)'GAS_CONSTANT' + write(10,900)rgas*1000/wm!real + + write(10,901)'TABLE_1'!See NOTE 8 + write(10,*)nt,np!integer + write(10,901)'TABLE_2' + write(10,*)nt,np!integer + write(10,901)'TABLE_3' + write(10,*)nt,np!integer + write(10,901)'TABLE_4' + write(10,*)nt,np!integer + write(10,901)'TABLE_5' + write(10,*)nt,np!integer + write(10,901)'TABLE_6' + write(10,*)nt,np!integer + write(10,901)'TABLE_7' + write(10,*)nt,np!integer + write(10,901)'TABLE_8'! See NOTE 9 + write(10,*)nt,np!integer + write(10,901)'TABLE_9' + write(10,*)nt,np!integer + write(10,901)'SAT_TABLE' + write(10,*)npsat,4,9 !integer + write(10,901)'$$SUPER_TABLE' + write(10,*)9 !integer (number of superheat tables, nn = 9) + + ! Output format units + 900 format (5ES17.7E3) + 901 format (A) + end do + + end subroutine writeRGPHeader + + + !---------------------------------------------------------------------! + ! SUPERHEAT TABLE ! + ! ! + ! writeSuperTable : Writes all super table items ! + ! writeSupTableHdr : Adds info about the super table ! + ! writeTable(...) : Writes the given data to file ! + ! ! + !---------------------------------------------------------------------! + + subroutine writeSuperTable() + call writeTable1() + call writeTable2() + call writeTable3() + call writeTable4() + call writeTable5() + call writeTable6() + call writeTable7() + call writeTable8() + call writeTable9() + end subroutine writeSuperTable + + subroutine writeSupTableHdr(tableno) + USE RGP_DATA + implicit none + character(len=1), intent(in) :: tableno + + write(10,901)'$TABLE_'//tableno + write(10,*)nt,np!integer(size of superheat arrays, See NOTE 10) + write(10,900) (t(i),i=1,nt) !|all following real + write(10,900) (p(j),j=1,np) !| + + ! Output format units + 900 format (5ES17.7E3) + 901 format (A) + end subroutine writeSupTableHdr + + subroutine writeTable1() + USE RGP_DATA + implicit none + + call writeSupTableHdr('1') + write(10,900) ((h(i,j),i=1,nt),j=1,np) !| + write(10,900) (tsat(i),i=1,nsat) !| + write(10,900) (hsat(i),i=1,nsat) !| + + ! Output format units + 900 format (5ES17.7E3) + end subroutine writeTable1 + + subroutine writeTable2() + USE RGP_DATA + implicit none + + call writeSupTableHdr('2') + write(10,900) ((w(i,j),i=1,nt),j=1,np) !| + write(10,900) (tsat(i),i=1,nsat) !| + write(10,900) (wsat(i),i=1,nsat) !| + + ! Output format units + 900 format (5ES17.7E3) + end subroutine writeTable2 + + subroutine writeTable3() + USE RGP_DATA + implicit none + + call writeSupTableHdr('3') + write(10,900) ((v(i,j),i=1,nt),j=1,np) !| + write(10,900) (tsat(i),i=1,nsat) !| + write(10,900) (vsat(i),i=1,nsat) !| + + ! Output format units + 900 format (5ES17.7E3) + end subroutine writeTable3 + + subroutine writeTable4() + USE RGP_DATA + implicit none + + call writeSupTableHdr('4') + write(10,900) ((cv(i,j),i=1,nt),j=1,np) !| + write(10,900) (tsat(i),i=1,nsat) !| + write(10,900) (cvsat(i),i=1,nsat) !| + + ! Output format units + 900 format (5ES17.7E3) + end subroutine writeTable4 + + subroutine writeTable5() + USE RGP_DATA + implicit none + + call writeSupTableHdr('5') + write(10,900) ((cp(i,j),i=1,nt),j=1,np) !| + write(10,900) (tsat(i),i=1,nsat) !| + write(10,900) (cpsat(i),i=1,nsat) !| + + ! Output format units + 900 format (5ES17.7E3) + end subroutine writeTable5 + + subroutine writeTable6() + USE RGP_DATA + implicit none + + call writeSupTableHdr('6') + write(10,900) ((dPdv(i,j),i=1,nt),j=1,np) !| + write(10,900) (tsat(i),i=1,nsat) !| + write(10,900) (dPdvsat(i),i=1,nsat) !| + + ! Output format units + 900 format (5ES17.7E3) + end subroutine writeTable6 + + subroutine writeTable7() + USE RGP_DATA + implicit none + + call writeSupTableHdr('7') + write(10,900) ((s(i,j),i=1,nt),j=1,np) !| + write(10,900) (tsat(i),i=1,nsat) !| + write(10,900) (ssat(i),i=1,nsat) !| + + ! Output format units + 900 format (5ES17.7E3) + end subroutine writeTable7 + + subroutine writeTable8() + USE RGP_DATA + implicit none + + call writeSupTableHdr('8') + write(10,900) ((eta(i,j),i=1,nt),j=1,np) !| + write(10,900) (tsat(i),i=1,nsat) !| + write(10,900) (etasat(i),i=1,nsat) !| + + ! Output format units + 900 format (5ES17.7E3) + end subroutine writeTable8 + + subroutine writeTable9() + USE RGP_DATA + implicit none + + call writeSupTableHdr('9') + write(10,900) ((tcx(i,j),i=1,nt),j=1,np) !| + write(10,900) (tsat(i),i=1,nsat) !| + write(10,900) (tcxsat(i),i=1,nsat) !| + + ! Output format units + 900 format (5ES17.7E3) + end subroutine writeTable9 + + !---------------------------------------------------------------------! + ! SATURATION TABLE ! + ! Note: 12.7 last paragraph: ! + ! The $$SAT_TABLE section does not have to exist under a $$$Database ! + ! access key, that the SUPERCOOLING parameter can be zero and that ! + ! only the $$SUPER_TABLE section is necessary. ! + ! ! + ! writeSatTable : Adds the sat table to the RGP file ! + ! writeSatData(...) : Writes the given data to file ! + ! ! + !---------------------------------------------------------------------! + + subroutine writeSatTable() + USE RGP_DATA + implicit none + + write (*,*) 'Sat Table' + write(10,901)'$$SAT_TABLE' + write(10,*)npsat,4,9 !(size of saturation arrays, see NOTE 11) + write(10,900) (psat2(i)*1000,i=1,npsat) + write(10,900) (tsat2(i),i=1,npsat) + write(10,900) (blank1(i),i=1,npsat) + write(10,900) (blank2(i),i=1,npsat) + + ! Liquid group + call writeSatData(hsatl, cpsatl, rhosatl, dPdrhosatl, ssatl, cvsatl,& + & wsatl, etasatl, tcxsatl) + + ! Vapour group + call writeSatData(hsatv, cpsatv, rhosatv, dPdrhosatv, ssatv, cvsatv, & + & wsatv, etasatv, tcxsatv) + + ! Output format units + 900 format (5ES17.7E3) + 901 format (A) + end subroutine writeSatTable + + subroutine writeSatData(H, Cp, Rho, dPdRho_T, S, Cv, SpSnd, eta, tcx) + implicit none + double precision, allocatable :: H(:), Cp(:), Rho(:), dPdRho_T(:), S(:), Cv(:), SpSnd(:), eta(:), tcx(:) + + write(10,900) (H(i),i=1,npsat) ! enthalpy + write(10,900) (Cp(i),i=1,npsat) ! Cp + write(10,900) (Rho(i),i=1,npsat) ! density + write(10,900) (dPdRho_T(i),i=1,npsat) ! change density with pressure const T + write(10,900) (S(i),i=1,npsat) ! entropy + write(10,900) (Cv(i),i=1,npsat) ! Cv + write(10,900) (SpSnd(i),i=1,npsat) ! Speed of Sound + write(10,900) (eta(i),i=1,npsat) ! molar viscosity + write(10,900) (tcx(i),i=1,npsat) ! thermal conductivity + + ! Output format units + 900 format (5ES17.7E3) + end subroutine writeSatData + + + !---------------------------------------------------------------------! + ! HELPER FUNCTIONS ! + !---------------------------------------------------------------------! + subroutine askQuestion(question, response) + implicit none + character(len=*), intent(in) :: question + double precision, intent(inout) :: response + + if (response .eq. 0) then + write(*,*) 'Enter - ', question + read(*,*) response + else + write(*,*) question, response + end if + end subroutine askQuestion + + subroutine askQuestion_int(question, response) + implicit none + character(len=*), intent(in) :: question + integer, intent(inout) :: response + + if (response .eq. 0) then + write(*,*) 'Enter - ', question + read(*,*) response + else + write(*,*) question, response + end if + end subroutine askQuestion_int + + subroutine getDblArg(position, value) + implicit none + character(len=255) :: arg + integer, intent(in) :: position + double precision, intent(out) :: value + + call get_command_argument(position, arg) + read(arg, *)value + end subroutine getDblArg + + subroutine getIntArg(position, value) + implicit none + character(len=255) :: arg + integer, intent(in) :: position + integer, intent(out) :: value + + call get_command_argument(position, arg) + read(arg, *)value + end subroutine getIntArg + + subroutine checkRPError(ierr, herr) + implicit none; + integer, intent(in) :: ierr + character(len=255), intent(in) :: herr + if (ierr .lt. 0) then + write (*,*) herr + call exit(-1) + end if + end subroutine checkRPError + + subroutine allocateArrays() + USE RGP_DATA + implicit none + + ! allocate 2D array + write (*,*) 'Allocating superheat array...' + ALLOCATE (t(nt), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 001' + ALLOCATE (p(np), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 002' + ALLOCATE (h(nt,np), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 003' + ALLOCATE (w(nt,np), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 004' + ALLOCATE (v(nt,np), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 005' + ALLOCATE (cv(nt,np), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 006' + ALLOCATE (cp(nt,np), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 007' + ALLOCATE (dPdv(nt,np), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 008' + ALLOCATE (s(nt,np), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 009' + ALLOCATE (eta(nt,np), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 010' + ALLOCATE (tcx(nt,np), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 011' + + ! allocate 1D array + nsat = np + ALLOCATE (tsat(nsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 012' + ALLOCATE (hsat(nsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 013' + ALLOCATE (wsat(nsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 014' + ALLOCATE (vsat(nsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 015' + ALLOCATE (cvsat(nsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 016' + ALLOCATE (cpsat(nsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 017' + ALLOCATE (dPdvsat(nsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 018' + ALLOCATE (ssat(nsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 019' + ALLOCATE (etasat(nsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 020' + ALLOCATE (tcxsat(nsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 021' + + ALLOCATE (tsat2(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 022' + ALLOCATE (psat2(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 023' + ALLOCATE (blank1(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 024' + ALLOCATE (blank2(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 025' + + ALLOCATE (hsatl(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 026' + ALLOCATE (wsatl(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 027' + ALLOCATE (rhosatl(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 028' + ALLOCATE (cvsatl(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 029' + ALLOCATE (cpsatl(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 030' + ALLOCATE (dPdrhosatl(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 031' + ALLOCATE (ssatl(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 032' + ALLOCATE (etasatl(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 033' + ALLOCATE (tcxsatl(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 034' + + ALLOCATE (hsatv(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 035' + ALLOCATE (wsatv(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 036' + ALLOCATE (rhosatv(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 037' + ALLOCATE (cvsatv(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 038' + ALLOCATE (cpsatv(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 039' + ALLOCATE (dPdrhosatv(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 040' + ALLOCATE (ssatv(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 041' + ALLOCATE (etasatv(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 042' + ALLOCATE (tcxsatv(npsat), STAT = AllocStatus) + if (AllocStatus /= 0) write (*,*)'out of memory 043' + end subroutine allocateArrays + +end program RGP diff --git a/RGPgenerator.F90 b/RGPgenerator.F90 deleted file mode 100644 index 2f74b97..0000000 --- a/RGPgenerator.F90 +++ /dev/null @@ -1,670 +0,0 @@ -! This program was written to generate Real Gas Property (RGP) files to be used -! with ANSYS CFX (Current version 14.5). The user is required to set certain -! input variables which generates an RGP file that follows the ANSYS guidelines -! see http://www.kxcad.net/ansys/ANSYS_CFX/help/help/Modelling/bk05ch10s07s01.html - -!Current Features: -! Works for CO2 -! Output file to ANSYS CFX version 14.5 specifications - -PROGRAM RGP - - !implicit double precision (a-h,o-z) - !implicit integer (i-k,m,n) - !implicit character*30 (l) - implicit NONE - integer ncmax - parameter (ncmax=20) !max number of components in mixture - double precision, dimension (ncmax) :: x(ncmax),xliq(ncmax),xvap(ncmax) - double precision, allocatable :: t(:),p(:),h(:,:),w(:,:),v(:,:),cv(:,:),cp(:,:),s(:,:),dPdv(:,:),eta(:,:),tcx(:,:), & - &tsat(:),hsat(:),wsat(:),vsat(:),cvsat(:),cpsat(:),ssat(:),dPdvsat(:),etasat(:),tcxsat(:), & - &tsat2(:),hsatl(:),wsatl(:),rhosatl(:),cvsatl(:),cpsatl(:),ssatl(:),dPdrhosatl(:),etasatl(:),tcxsatl(:), & - &psat2(:),hsatv(:),wsatv(:),rhosatv(:),cvsatv(:),cpsatv(:),ssatv(:),dPdrhosatv(:),etasatv(:),tcxsatv(:), & - &blank1(:),blank2(:),arr(:) - double precision :: acf,conp,conv,conpl,convl,conpv,convv,cpcrit,cvcrit,density,dpdrho,q,rgas,rhocrit,rhol,rhov, & - &dc,dip,dl,dm,dmax,dpdrhol,dpdrhov,dv,e,enthalpy,enthalpyl,enthalpyv,pstep,pstepsat,scrit,speedcrit, & - &entropy,entropyl,entropyv,hcrit,htyp,pc,pm,pmax,pmaxsat,pminsat,psat,preslower,upperPres,speedOsound, & - &speedOsoundl,speedOsoundv,supercool,tc,templower,uppertemp,TminSat,TmaxSat,thermalconductivity,thermalconductivityl, & - &thermalconductivityv,thermcrit,tliq,tm,tmax,tmin,tnbp,tpp,ttp,tstep,tvap,visccrit,viscosity,viscosityl,viscosityv,wm, & - &y,zc,dbl - integer :: i,j,ierr,kq,np,nt,nsat,npsat,ilng - character :: hrf*3, herr*255, lfluid*30, loutfile*30 - character :: lFileName*255, hnam*12, hn80*80, hstr*255 - character :: hcasn*12 - character(LEN=255) :: hf(ncmax) - character*7 hfmix - integer :: prevStatus - integer :: currStatus - integer :: tenPercent - integer :: AllocStatus - i=1 - !PROMPT FOR OUTPUT FILE - WRITE(*,*) ' Carleton University' - WRITE(*,*) ' Department of Mechanical and Aerospace Engineering' - WRITE(*,*) - WRITE(*,*) ' 88888888ba ,ad8888ba, 88888888ba ' - WRITE(*,*) ' 88 "8b d8" "8b 88 "8b ' - WRITE(*,*) ' 88 .8P d8 88 8P ' - WRITE(*,*) ' 88aaaaaa8P 88 88aaaaaa8P ' - WRITE(*,*) ' 88""""88 88 88888 88"""""" ' - WRITE(*,*) ' 88 `8b Y8, `88 88 ' - WRITE(*,*) ' 88 `8b Y8a. .a88 88 ' - WRITE(*,*) ' 88 `8b `"Y88888P" 88 ' - WRITE(*,*) - WRITE(*,*) ' FILE GENERATOR' - WRITE(*,*) - WRITE(*,'(A)') 'Output Data File Name (i.e. filename.txt):' - READ (*,*) lOutfile - OPEN(UNIT=10,FILE=lOutFile,STATUS='NEW') - lfluid='SCO2' - WRITE(*,'(A)') 'Fluid file name (i.e. filename.fld):' - READ (*,*) lFileName - hf(1)=lFileName - hfmix='hmx.bnc' - hrf='DEF' - CALL SETUP (i,hf,hfmix,hrf,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - CALL INFO (1,wm,ttp,tnbp,tc,pc,dc,zc,acf,dip,rgas) - CALL NAME (1,hnam,hn80,hcasn) - Call PASSCMN ('ptpn',0,1,0, hstr,ilng,tpp,arr,ierr,herr) - dc=wm*dc!Convert density to kg/m^3 - Dmax=wm*Dmax - !tpp = 517.95 !kPa *********** ADDED FROM CO2.FLD ******* - WRITE(*,*) - WRITE(*,*) 'hnam',hnam - WRITE(*,*) 'hn80',hn80 - WRITE(*,*) 'Molecular Weight: ',wm,'g/mol' -! WRITE(*,*) 'Accentric Factor [-]: ',acf SETUP error -! WRITE(*,*) 'Dipole Moment [debye]: ',dip - WRITE(*,*) 'Triple Point Temperature: ',ttp,'K' - WRITE(*,*) 'Triple Point Pressure: ',tpp,'kPa' - WRITE(*,*) 'Normal Boiling Point: ',tnbp,'K' - - WRITE(*,*) 'Critical Temperature: ',tc,'K' - WRITE(*,*) 'Critical Pressure: ',pc,'kPa' - WRITE(*,*) 'Critical Density: ',dc,'kg/m^3' - WRITE(*,*) 'Compressibility at Critical Point: ',zc - WRITE(*,*) 'Gas Constant: ',rgas,'J/mol-K' - CALL MAXT (x,tm,pm,Dm,ierr,herr) - CALL LIMITS(htyp,x,tmin,tmax,Dmax,pmax) - WRITE(*,*) - WRITE(*,*)'Limits for this fluid:' - WRITE(*,*) - WRITE(*,*)'Min Temperature: ',tmin,'K' - WRITE(*,*)'Max Temperature: ',tmax,'K' - WRITE(*,*)'Max Saturation temperature: ',tm,'K' - WRITE(*,*)'Max Density: ',Dmax,'kg/m^3' - WRITE(*,*)'Max Pressure: ',pmax,'kPa' - - ! PROMPT FOR BOUNDS AND STEP SIZE - - WRITE(*,*)'--------------------------------------------' - WRITE(*,*) - WRITE(*,*)'Enter - Lower bound for Temperature [K]:' - READ(*,*)TempLower - WRITE(*,*)'Enter - Upper bound for Temperature [K]:' - READ(*,*)UpperTemp - WRITE(*,*)'Enter - Number of Steps for Temperature:' - READ(*,*)nt - WRITE(*,*)'Enter - Lower bound for Pressure [Pa]:' - READ(*,*)PresLower - WRITE(*,*)'Enter - Upper bound for Pressure [Pa]:' - READ(*,*)UpperPres - WRITE(*,*)'Enter - Number of Steps for Pressure:' - READ(*,*)np - WRITE(*,*)'Enter - Minimum Temperature for Saturation Tables:' - READ(*,*)TminSat - WRITE(*,*)'Enter - Maximum Temperature for Saturation Tables:' - READ(*,*)TmaxSat - WRITE(*,*)'Enter - Number of Steps for Saturation Tables:' - READ(*,*)npsat - - SuperCool=0 - - WRITE(*,*) - WRITE(*,*)'Generating RGP File Please Wait...' - WRITE(*,*) - - ! SETUP OUTPUT FILE - WRITE(10,901)' $$$$HEADER' - WRITE(10,901)' $$$',hnam!,lfluid !character*8 (key into the RGP file see NOTE 3) if this program is used with a fluid other than CO2, this line will need to be changed, along with the corresponding line in the $$$DATA header - WRITE(10,*) 1!enter an integer (this line is ignored in ANSYS CFX) - WRITE(10,901)'$$PARAM'! See NOTE 4 - WRITE(10,*) 26!enter an integer(number of parameters) - WRITE(10,901)'DESCRIPTION' - WRITE(10,901)'Carbon Dioxide from NIST'!character*50(description of the material) - WRITE(10,901)'NAME' - WRITE(10,901)hnam!lfluid!character*8(material name, same as $$$) - WRITE(10,901)'INDEX' - WRITE(10,901)'SCO2'!lfluid!character*50 (index into clients RGDB program) - WRITE(10,901)'DATABASE' - WRITE(10,901)'NIST REAL GAS PROPERTY DATABASE' - WRITE(10,901)'MODEL' - WRITE(10,*)3!integer(level of property info available 1,2,or 3) - WRITE(10,901)'UNITS' - WRITE(10,*)1!integer(unit system of 1 ,2,3,4 or 5) - WRITE(10,901)'PMIN_SUPERHEAT'! at Triple Point Pressure? - WRITE(10,900)PresLower!real (Pa) - WRITE(10,901)'PMAX_SUPERHEAT'! at Critical Pressure? - WRITE(10,900)UpperPres!real (Pa) - WRITE(10,901)'TMIN_SUPERHEAT'! at Triple Point Temperature? - WRITE(10,900)TempLower!real (K) - WRITE(10,901)'TMAX_SUPERHEAT'! at Critical Temperature? - WRITE(10,900)UpperTemp!real (K) - WRITE(10,901)'SUPERCOOLING' - WRITE(10,900)SuperCool!real(supercooling level in superheat tables)(Optional) - WRITE(10,901)'P_CRITICAL' - WRITE(10,900)Pc*1000!real - WRITE(10,901)'P_TRIPLE'! See NOTE 7 - WRITE(10,900)tpp*1000!real - WRITE(10,901)'T_CRITICAL' - WRITE(10,900)Tc!real - WRITE(10,901)'T_TRIPLE' - WRITE(10,900)Ttp!real - WRITE(10,901)'GAS_CONSTANT' - WRITE(10,900)rgas*1000/wm!real - - WRITE(10,901)'TABLE_1'!See NOTE 8 - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_2' - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_3' - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_4' - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_5' - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_6' - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_7' - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_8'! See NOTE 9 - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_9' - WRITE(10,*)nt,np!integer - WRITE(10,901)'SAT_TABLE' - WRITE(10,*)npsat,4,9 !integer - - WRITE(10,901) '$$$$DATA' - WRITE(10,901) '$$$SCO2'!,lfluid !character*8 (key into the RGP file see NOTE 3) - WRITE(10,*) 1!enter an integer (this line is ignored in ANSYS CFX) - WRITE(10,901)'$$PARAM'! See NOTE 4 - WRITE(10,*) 26!enter an integer(number of parameters) - WRITE(10,901)'DESCRIPTION' - WRITE(10,901)'Carbon Dioxide from NIST'!character*50(description of the material) - WRITE(10,901)'NAME' - WRITE(10,901)'SCO2'!lfluid!character*8(material name, same as $$$) - WRITE(10,901)'INDEX' - WRITE(10,901)'SCO2'!lfluid!character*50 (index into clients RGDB program) - WRITE(10,901)'DATABASE' - WRITE(10,901)'NIST REAL GAS PROPERTY DATABASE' - WRITE(10,901)'MODEL' - WRITE(10,*)3!integer(level of property info available 1,2,or 3) - WRITE(10,901)'UNITS' - WRITE(10,*)1!integer(unit system of 1 ,2,3,4 or 5) - WRITE(10,901)'PMIN_SUPERHEAT'! at Triple Point Pressure? - WRITE(10,900)PresLower!real (Pa) - WRITE(10,901)'PMAX_SUPERHEAT'! at Critical Pressure? - WRITE(10,900)UpperPres!real (Pa) - WRITE(10,901)'TMIN_SUPERHEAT'! at Triple Point Temperature? - WRITE(10,900)TempLower!real (K) - WRITE(10,901)'TMAX_SUPERHEAT'! at Critical Temperature? - WRITE(10,900)UpperTemp!real (K) - WRITE(10,901)'SUPERCOOLING' - WRITE(10,900)SuperCool!real(supercooling level in superheat tables)(Optional) - WRITE(10,901)'P_CRITICAL' - WRITE(10,900)Pc*1000!real - WRITE(10,901)'P_TRIPLE'! See NOTE 7 - WRITE(10,900)tpp*1000!real - WRITE(10,901)'T_CRITICAL' - WRITE(10,900)Tc!real - WRITE(10,901)'T_TRIPLE' - WRITE(10,900)Ttp!real - WRITE(10,901)'GAS_CONSTANT' - WRITE(10,900)rgas*1000/wm!real - - WRITE(10,901)'TABLE_1'!See NOTE 8 - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_2' - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_3' - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_4' - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_5' - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_6' - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_7' - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_8'! See NOTE 9 - WRITE(10,*)nt,np!integer - WRITE(10,901)'TABLE_9' - WRITE(10,*)nt,np!integer - WRITE(10,901)'SAT_TABLE' - WRITE(10,*)npsat,4,9 !integer - - WRITE(10,901)'$$SUPER_TABLE' - WRITE(10,*)9 !integer (number of superheat tables, nn = 9) - tstep = (UpperTemp-TempLower)/REAL(nt-1) - pstep = (UpperPres-PresLower)/REAL(np-1) - - ! allocate 2D array - WRITE (*,*) 'Allocating superheat array...' - ALLOCATE (t(nt), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 001' - ALLOCATE (p(np), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 002' - ALLOCATE (h(nt,np), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 003' - ALLOCATE (w(nt,np), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 004' - ALLOCATE (v(nt,np), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 005' - ALLOCATE (cv(nt,np), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 006' - ALLOCATE (cp(nt,np), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 007' - ALLOCATE (dPdv(nt,np), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 008' - ALLOCATE (s(nt,np), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 009' - ALLOCATE (eta(nt,np), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 010' - ALLOCATE (tcx(nt,np), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 011' - - ! allocate 1D array - nsat = np - ALLOCATE (tsat(nsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 012' - ALLOCATE (hsat(nsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 013' - ALLOCATE (wsat(nsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 014' - ALLOCATE (vsat(nsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 015' - ALLOCATE (cvsat(nsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 016' - ALLOCATE (cpsat(nsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 017' - ALLOCATE (dPdvsat(nsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 018' - ALLOCATE (ssat(nsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 019' - ALLOCATE (etasat(nsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 020' - ALLOCATE (tcxsat(nsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 021' - - CALL TPFLSH (tc,pc,x,rhoCrit,dl,dv,xliq,xvap,q,e,Hcrit,Scrit,cvCrit,cpCrit,SpeedCrit,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - CALL DPDD (tc,rhoCrit,x,dPdrhov) - IF (ierr.ne.0) write (*,*) herr - CALL TRNPRP (tc,rhoCrit,x,viscCrit,ThermCrit,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - - !Fill superheat arrays with RGP data - !For 3D table: - !any point in the liquid area gets filled with saturation properties at that pressure - !any point below the crit temp but above the crit pressure gets filled with phi(T, P) (ie normally) - !For 2D tables: - !for every pressure in the 2D table there is an entry in the 1D table - !for pressures above Pcrit, the entry is of phi(Tmin, P) instead of phi(Tsat, P) - write (*,*) 'Filling superheat table...' - prevStatus = 0 - currStatus = 0 - tenPercent = INT(UpperPres-PresLower)/10 - j=1 - p(j) = PresLower - DO WHILE (j <= np) - i=1 - t(i) = TempLower - - ! show progress every 10% - currStatus = INT(p(j) - PresLower)/tenPercent - IF (prevStatus /= currStatus) THEN - write (*,*) (currStatus*10), '%' - prevStatus = currStatus - END IF - - DO WHILE (i <= nt) - IF (p(j) <= pc*1000) THEN - ! find the sat temp for this pressure, the if we are above it, work normal - ! if we are below the sat temp, then fill 2D table with phi(Tsat, p) - PSat= p(j)/1000 - q=1 - CALL PQFLSH (PSat,q,x,kq,tvap,rhov,Dl,Dv,xliq,xvap,e,enthalpyv,entropyv,convv,conpv,speedOsoundv,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - - IF (t(i) <= tvap) THEN - ! use Tsat properties at the given p - h(i,j) = enthalpyv*1000/wm - w(i,j) = speedOsoundv - v(i,j) = 1/(rhov*wm) - cv(i,j) = convv*1000/wm - cp(i,j) = conpv*1000/wm - s(i,j) = entropyv*1000/wm - CALL DPDD (tvap,rhov,x,dPdrhov) - IF (ierr.ne.0) write (*,*) herr - dPdv(i,j) = (dPdrhov/wm*1000*(-(rhov*wm)**2)) - CALL TRNPRP (tvap,rhov,x,viscosityv,thermalconductivityv,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - eta(i,j) = viscosityv/1000000 - tcx(i,j) = thermalconductivityv - ELSE - ! get properties normally - CALL TPFLSH (t(i),(p(j)/1000),x,density,dl,dv,xliq,xvap,q,e,enthalpy,entropy,conv,conp,speedOsound,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - h(i,j) = enthalpy*1000/wm - w(i,j) = speedOsound - v(i,j) = 1/(density*wm) - cv(i,j) = conv*1000/wm - cp(i,j) = conp*1000/wm - s(i,j) = entropy*1000/wm - CALL DPDD (t(i),density,x,dPdrho) - IF (ierr.ne.0) write (*,*) herr - dPdv(i,j) = (dPdrho/wm*1000*(-(density*wm)**2)) - CALL TRNPRP (t(i),density,x,viscosity,thermalconductivity,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - eta(i,j) = viscosity/1000000 - tcx(i,j) = thermalconductivity - END IF - ELSE - ! P > Pcrit - - ! run normally - CALL TPFLSH (t(i),(p(j)/1000),x,density,dl,dv,xliq,xvap,q,e,enthalpy,entropy,conv,conp,speedOsound,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - h(i,j) = enthalpy*1000/wm - w(i,j) = speedOsound - v(i,j) = 1/(density*wm) - cv(i,j) = conv*1000/wm - cp(i,j) = conp*1000/wm - s(i,j) = entropy*1000/wm - CALL DPDD (t(i),density,x,dPdrho) - IF (ierr.ne.0) write (*,*) herr - dPdv(i,j) = (dPdrho/wm*1000*(-(density*wm)**2)) - CALL TRNPRP (t(i),density,x,viscosity,thermalconductivity,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - eta(i,j) = viscosity/1000000 - tcx(i,j) = thermalconductivity - END IF - - i = i + 1 - IF (i <= nt) THEN - t(i) = t(i-1)+tstep - END IF - END DO - - ! Fill 1D Table (Sat line) - IF (p(j)<(Pc*1000))THEN - ! fill with phi(Tsat, P) - PSat= p(j)/1000 - q=1 - CALL PQFLSH (PSat,q,x,kq,tvap,rhov,Dl,Dv,xliq,xvap,e,enthalpyv,entropyv,convv,conpv,speedOsoundv,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - tsat(j) = tvap - hsat(j) = enthalpyv*1000/wm - wsat(j) = speedOsoundv - vsat(j) = 1/(rhov*wm) - cvsat(j) = convv*1000/wm - cpsat(j) = conpv*1000/wm - ssat(j) = entropyv*1000/wm - CALL DPDD (tvap,rhov,x,dPdrhov) - IF (ierr.ne.0) write (*,*) herr - dPdvsat(j) = (dPdrhov/wm*1000*(-(rhov*wm)**2)) - CALL TRNPRP (tvap,rhov,x,viscosityv,thermalconductivityv,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - etasat(j) = viscosityv/1000000 - tcxsat(j) = thermalconductivityv - ELSE - ! fill with phi(Tmin, P) - PSat= p(j)/1000 - q=1 - CALL TPFLSH (TempLower,(p(j)/1000),x,density,dl,dv,xliq,xvap,q,e,enthalpy,entropy,conv,conp,speedOsound,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - tsat(j) = TempLower - hsat(j) = enthalpy*1000/wm - wsat(j) = speedOsound - vsat(j) = 1/(density*wm) - cvsat(j) = conv*1000/wm - cpsat(j) = conp*1000/wm - ssat(j) = entropy*1000/wm - CALL DPDD (TempLower,density,x,dPdrho) - IF (ierr.ne.0) write (*,*) herr - dPdvsat(j) = (dPdrho/wm*1000*(-(density*wm)**2)) - CALL TRNPRP (TempLower,density,x,viscosity,thermalconductivity,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - etasat(j) = viscosity/1000000 - tcxsat(j) = thermalconductivity - END IF - - j = j+1 - IF (j <= np) THEN - p(j) = p(j-1)+pstep - END IF - END DO - - q=1 - CALL TQFLSH (TminSat,q,x,kq,PminSat,density,Dl,Dv,x,y,e,enthalpy,entropy,conv,conp,speedOsound,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - CALL TQFLSH (TmaxSat,q,x,kq,PmaxSat,density,Dl,Dv,x,y,e,enthalpy,entropy,conv,conp,speedOsound,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - pstepsat = (Pmaxsat-Pminsat)/REAL(npsat-1) - - WRITE (*,*) 'Allocating saturation table...' - ALLOCATE (tsat2(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 022' - ALLOCATE (psat2(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 023' - ALLOCATE (blank1(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 024' - ALLOCATE (blank2(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 025' - - ALLOCATE (hsatl(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 026' - ALLOCATE (wsatl(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 027' - ALLOCATE (rhosatl(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 028' - ALLOCATE (cvsatl(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 029' - ALLOCATE (cpsatl(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 030' - ALLOCATE (dPdrhosatl(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 031' - ALLOCATE (ssatl(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 032' - ALLOCATE (etasatl(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 033' - ALLOCATE (tcxsatl(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 034' - - ALLOCATE (hsatv(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 035' - ALLOCATE (wsatv(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 036' - ALLOCATE (rhosatv(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 037' - ALLOCATE (cvsatv(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 038' - ALLOCATE (cpsatv(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 039' - ALLOCATE (dPdrhosatv(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 040' - ALLOCATE (ssatv(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 041' - ALLOCATE (etasatv(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 042' - ALLOCATE (tcxsatv(npsat), STAT = AllocStatus) - IF (AllocStatus /= 0) write (*,*)'out of memory 043' - - ! fill saturation table - WRITE (*,*) 'Filling saturation table...' - j=1 - psat2(j) = PminSat - DO WHILE (j <= npsat) - q=1 - CALL PQFLSH (psat2(j),q,x,kq,tvap,rhov,Dl,Dv,xliq,xvap,e,enthalpyv,entropyv,convv,conpv,speedOsoundv,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - tsat2(j) = tvap - blank1(j)=0 - blank2(j)=0 - hsatv(j) = enthalpyv*1000/wm - wsatv(j) = speedOsoundv - rhosatv(j) = (rhov*wm) - cvsatv(j) = convv*1000/wm - cpsatv(j) = conpv*1000/wm - ssatv(j) = entropyv*1000/wm - CALL DPDD (tvap,rhov,x,dPdrhov) - IF (ierr.ne.0) write (*,*) herr - dPdrhosatv(j) = wm/(dPdrhov*1000) - CALL TRNPRP (tsat2(j),rhov,x,viscosityv,thermalconductivityv,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - etasatv(j) = viscosityv/1000000 - tcxsatv(j) = thermalconductivityv - - q=0 - CALL PQFLSH (psat2(j),q,x,kq,tliq,rhol,Dl,Dv,xliq,xvap,e,enthalpyl,entropyl,convl,conpl,speedOsoundl,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - hsatl(j) = enthalpyl*1000/wm - wsatl(j) = speedOsoundl - rhosatl(j) = (rhol*wm) - cvsatl(j) = convl*1000/wm - cpsatl(j) = conpl*1000/wm - ssatl(j) = entropyl*1000/wm - CALL DPDD (tliq,rhol,x,dPdrhol) - IF (ierr.ne.0) write (*,*) herr - dPdrhosatl(j) = wm/(dPdrhol*1000) - CALL TRNPRP (tsat2(j),rhol,x,viscosityl,thermalconductivityl,ierr,herr) - IF (ierr.ne.0) write (*,*) herr - etasatl(j) = viscosityl/1000000 - tcxsatl(j) = thermalconductivityl - - j = j+1 - IF (j <= npsat) THEN - psat2(j) = psat2(j-1) + pstepsat - END IF - END DO - - WRITE (*,*) 'Writing data to file...' - !Print data to the text file in the correct format - WRITE (*,*) 'Table 1' - WRITE(10,901)'$TABLE_1' - WRITE(10,*)nt,np!integer(size of superheat arrays, See NOTE 10) - WRITE(10,900) (t(i),i=1,nt) !|all following real - WRITE(10,900) (p(j),j=1,np) !| - WRITE(10,900) ((h(i,j),i=1,nt),j=1,np) !| - WRITE(10,900) (tsat(i),i=1,nsat) !| - WRITE(10,900) (hsat(i),i=1,nsat) !| - - WRITE (*,*) 'Table 2' - WRITE(10,901)'$TABLE_2' - WRITE(10,*)nt,np!integer(size of superheat arrays, See NOTE 10) - WRITE(10,900) (t(i),i=1,nt) !|all following real - WRITE(10,900) (p(j),j=1,np) !| - WRITE(10,900) ((w(i,j),i=1,nt),j=1,np) !| - WRITE(10,900) (tsat(i),i=1,nsat) !| - WRITE(10,900) (wsat(i),i=1,nsat) !| - - WRITE (*,*) 'Table 3' - WRITE(10,901)'$TABLE_3' - WRITE(10,*)nt,np!integer(size of superheat arrays, See NOTE 10) - WRITE(10,900) (t(i),i=1,nt) !|all following real - WRITE(10,900) (p(j),j=1,np) !| - WRITE(10,900) ((v(i,j),i=1,nt),j=1,np) !| - WRITE(10,900) (tsat(i),i=1,nsat) !| - WRITE(10,900) (vsat(i),i=1,nsat) !| - - WRITE (*,*) 'Table 4' - WRITE(10,901)'$TABLE_4' - WRITE(10,*)nt,np!integer(size of superheat arrays, See NOTE 10) - WRITE(10,900) (t(i),i=1,nt) !|all following real - WRITE(10,900) (p(j),j=1,np) !| - WRITE(10,900) ((cv(i,j),i=1,nt),j=1,np) !| - WRITE(10,900) (tsat(i),i=1,nsat) !| - WRITE(10,900) (cvsat(i),i=1,nsat) !| - - WRITE (*,*) 'Table 5' - WRITE(10,901)'$TABLE_5' - WRITE(10,*)nt,np!integer(size of superheat arrays, See NOTE 10) - WRITE(10,900) (t(i),i=1,nt) !|all following real - WRITE(10,900) (p(j),j=1,np) !| - WRITE(10,900) ((cp(i,j),i=1,nt),j=1,np) !| - WRITE(10,900) (tsat(i),i=1,nsat) !| - WRITE(10,900) (cpsat(i),i=1,nsat) !| - - WRITE (*,*) 'Table 6' - WRITE(10,901)'$TABLE_6' - WRITE(10,*)nt,np!integer(size of superheat arrays, See NOTE 10) - WRITE(10,900) (t(i),i=1,nt) !|all following real - WRITE(10,900) (p(j),j=1,np) !| - WRITE(10,900) ((dPdv(i,j),i=1,nt),j=1,np) !| - WRITE(10,900) (tsat(i),i=1,nsat) !| - WRITE(10,900) (dPdvsat(i),i=1,nsat) !| - - WRITE (*,*) 'Table 7' - WRITE(10,901)'$TABLE_7' - WRITE(10,*)nt,np!integer(size of superheat arrays, See NOTE 10) - WRITE(10,900) (t(i),i=1,nt) !|all following real - WRITE(10,900) (p(j),j=1,np) !| - WRITE(10,900) ((s(i,j),i=1,nt),j=1,np) !| - WRITE(10,900) (tsat(i),i=1,nsat) !| - WRITE(10,900) (ssat(i),i=1,nsat) !| - - WRITE (*,*) 'Table 8' - WRITE(10,901)'$TABLE_8' - WRITE(10,*)nt,np!integer(size of superheat arrays, See NOTE 10) - WRITE(10,900) (t(i),i=1,nt) !|all following real - WRITE(10,900) (p(i),i=1,np) !| - WRITE(10,900) ((eta(i,j),i=1,nt),j=1,np) !| - WRITE(10,900) (tsat(i),i=1,nsat) !| - WRITE(10,900) (etasat(i),i=1,nsat) !| - - WRITE (*,*) 'Table 9' - WRITE(10,901)'$TABLE_9' - WRITE(10,*)nt,np!integer(size of superheat arrays, See NOTE 10) - WRITE(10,900) (t(i),i=1,nt) !|all following real - WRITE(10,900) (p(j),j=1,np) !| - WRITE(10,900) ((tcx(i,j),i=1,nt),j=1,np) !| - WRITE(10,900) (tsat(i),i=1,nsat) !| - WRITE(10,900) (tcxsat(i),i=1,nsat) !| - - WRITE (*,*) 'Sat Table' - WRITE(10,901)'$$SAT_TABLE' - WRITE(10,*)npsat,4,9 !(size of saturation arrays, see NOTE 11) - WRITE(10,900) (psat2(i)*1000,i=1,npsat) - WRITE(10,900) (tsat2(i),i=1,npsat) - WRITE(10,900) (blank1(i),i=1,npsat) - WRITE(10,900) (blank2(i),i=1,npsat) - WRITE(10,900) (hsatl(i),i=1,npsat) - WRITE(10,900) (cpsatl(i),i=1,npsat) - WRITE(10,900) (rhosatl(i),i=1,npsat) - WRITE(10,900) (dPdrhosatl(i),i=1,npsat)! - WRITE(10,900) (ssatl(i),i=1,npsat) ! - WRITE(10,900) (cvsatl(i),i=1,npsat)! - WRITE(10,900) (wsatl(i),i=1,npsat) - WRITE(10,900) (etasatl(i),i=1,npsat) - WRITE(10,900) (tcxsatl(i),i=1,npsat) - WRITE(10,900) (hsatv(i),i=1,npsat) - WRITE(10,900) (cpsatv(i),i=1,npsat)! - WRITE(10,900) (rhosatv(i),i=1,npsat) - WRITE(10,900) (dPdrhosatv(i),i=1,npsat)! - WRITE(10,900) (ssatv(i),i=1,npsat) - WRITE(10,900) (cvsatv(i),i=1,npsat) - WRITE(10,900) (wsatv(i),i=1,npsat) - WRITE(10,900) (etasatv(i),i=1,npsat) - WRITE(10,900) (tcxsatv(i),i=1,npsat) - - CLOSE(10) - - WRITE(*,*) - WRITE(*,*)'RGP File Generated!' - READ(*,*) - 900 format (5ES17.7E3) - 901 format (A) - - -END PROGRAM RGP diff --git a/fluids/copy fluids here.txt b/fluids/copy fluids here.txt new file mode 100644 index 0000000..e69de29 diff --git a/fortran/copy fortran here.txt b/fortran/copy fortran here.txt new file mode 100644 index 0000000..e69de29 diff --git a/mixtures/copy mixtures here.txt b/mixtures/copy mixtures here.txt new file mode 100644 index 0000000..e69de29 diff --git a/patchUtility.py b/patchUtility.py index d2cf162..aee47cf 100644 --- a/patchUtility.py +++ b/patchUtility.py @@ -1,3 +1,4 @@ +from __future__ import print_function import sys original = sys.argv[1] @@ -8,12 +9,17 @@ if (l) ptpn(icomp)=dbl elseif (hvar.eq.'txeos') then""" -with open(original, 'r') as utility_file: - utility = utility_file.read(); - -utility = utility.replace("elseif (hvar.eq.'txeos') then",newcode) +def addPTPN(): + with open(original, 'r') as utility_file: + utility = utility_file.read() + + # don't patch if already patched + # (unlikely as we get the original file) + if "elseif (hvar.eq.'ptpn') then" not in utility: + utility = utility.replace("elseif (hvar.eq.'txeos') then",newcode) -with open(dest, 'w') as utility_file: - utility_file.write(utility) + with open(dest, 'w') as utility_file: + utility_file.write(utility) -print('patch complete') +addPTPN() +print('REFPROP: Patch complete', end="") \ No newline at end of file diff --git a/test/check_format.py b/test/check_format.py new file mode 100644 index 0000000..7bfef7e --- /dev/null +++ b/test/check_format.py @@ -0,0 +1,44 @@ +import os +from math import floor + +def loadRGP(rgp_name): + with open(rgp_name, "r") as rgp_fp: + rgp_data = rgp_fp.read() + + # parse example + rgp_data = rgp_data.split('$$$$') + data_table = rgp_data[2].split('$TABLE_') + sat_table = data_table[9].split('$$SAT_TABLE') + data_table[9] = sat_table[0] + sat_table = sat_table[1] + data_table = [data.replace('\n', '').split() for data in data_table[1:len(data_table)]] + + + for i in range(0, len(data_table)): + data = data_table[i][3:len(data_table[i])] + nPoints = int(data_table[i][2]) + data_table[i] = [] + for j in range(0, floor(len(data)/nPoints)): + data_table[i].append(data[j*nPoints:(j+1)*nPoints]) + + return data_table + +example_rgp = loadRGP('format_example.rgp') +our_rgp = loadRGP('out.rgp') + +for i in range(1, len(example_rgp)): + table_ex = example_rgp[i] + table_ou = our_rgp[i] + + for j in range(1, len(table_ex)): + phi_ex = [float(num) for num in table_ex[j]] + phi_ou = [float(num) for num in table_ou[j]] + + for k in range(1, len(phi_ex)): + diff = (phi_ex[k] - phi_ou[k]) / phi_ex[k] + + if diff>1e-2: + print('Table', i+1, 'phi', j+1, 'item', k+1, 'diff = ', round(diff,4), + 'ex', round(phi_ex[k],4), 'ou', round(phi_ou[k],4)) + +print('hi') \ No newline at end of file diff --git a/test/format_example.rgp b/test/format_example.rgp new file mode 100644 index 0000000..473303c --- /dev/null +++ b/test/format_example.rgp @@ -0,0 +1,282 @@ +$$$$HEADER +$$$R134A +1 +$$PARAM +26 +DESCRIPTION +Refrigerant R134a from NIST +NAME +R134A +INDEX +R134A +MODEL +2 +UNITS +1 +PMIN_SUPERHEAT +55000.00000000000 +PMAX_SUPERHEAT +400000.0000000000 +TMIN_SUPERHEAT +240.0000000000000 +TMAX_SUPERHEAT +400.0000000000000 +TMIN_SATURATION +240.0000000000000 +TMAX_SATURATION +350.0000000000000 +SUPERCOOLING +0.0000000000000000E+00 +P_CRITICAL +405928.0E+3 +P_TRIPLE +0.0000000 +T_CRITICAL +374.21 +T_TRIPLE +169.85 +TABLE_1 +5 5 +TABLE_2 +5 5 +TABLE_3 +5 5 +TABLE_4 +5 5 +TABLE_5 +5 5 +TABLE_6 +5 5 +TABLE_7 +5 5 +TABLE_8 +5 5 +TABLE_9 +5 5 +SAT_TABLE +5 4 9 +$$$$DATA +$$$R134A +1 +$$PARAM +26 +DESCRIPTION + Refrigerant R134a from NIST +NAME +R134A +INDEX +R134A +MODEL +2 +UNITS +1 +PMIN_SUPERHEAT +55000.00000000000 +PMAX_SUPERHEAT +400000.0000000000 +TMIN_SUPERHEAT +240.0000000000000 +TMAX_SUPERHEAT +400.0000000000000 +TMIN_SATURATION +240.0000000000000 +TMAX_SATURATION +350.0000000000000 +SUPERCOOLING +0.0000000000000000E+00 +P_CRITICAL +405928.0E+3 +P_TRIPLE +0.00000000 +T_CRITICAL +374.21 +T_TRIPLE +169.85 +TABLE_1 +5 5 +TABLE_2 +5 5 +TABLE_3 +5 5 +TABLE_4 +5 5 +TABLE_5 +5 5 +TABLE_6 +5 5 +TABLE_7 +5 5 +TABLE_8 +5 5 +TABLE_9 +5 5 +SAT_TABLE +5 4 9 +$$SUPER_TABLE +9 +$TABLE_1 +5 5 + 240.0000000000000 280.0000000000000 320.0000000000000 360.0000000000000 400.0000000000000 + 55000.00000000000 141250.0000000000 227500.0000000000 313750.0000000000 400000.0000000000 + 380537.4694438053 411338.2663531548 444977.3409773413 481319.1776109927 520216.5694145825 + 389653.7412125674 409809.2455431806 443799.3654908382 480400.4939131565 519497.1970340972 + 397134.4582379403 408220.8885809304 442595.4896125616 479470.0507244235 518772.4705830857 + 402478.8117714382 406565.3815292667 441363.8343756341 478527.3570462614 518042.2612321422 + 406650.5418732642 406650.5418732642 440102.2828142202 477571.8870479560 517306.4350445317 + 234.5147959861216 254.6474401194407 266.4385386783982 275.1224493134379 282.1192618430768 + 376541.6333194062 389653.7412125674 397134.4582379403 402478.8117714382 406650.5418732642 +$TABLE_2 + 5 5 + 240.0000000000000 280.0000000000000 320.0000000000000 360.0000000000000 400.0000000000000 + 55000.00000000000 141250.0000000000 227500.0000000000 313750.0000000000 400000.0000000000 + 146.0231331684759 157.8078814382528 168.5506303571756 178.5170722927099 187.8737728104760 + 147.0869026517344 155.2427165817653 166.8523505651282 177.3642207508735 187.0763724397798 + 147.8329014568262 152.5632728477825 165.1148711540457 176.1977255190386 186.2743889887011 + 147.9104560214784 149.7535601390501 163.3352394640403 175.0170180126285 185.4677211216712 + 147.6607838596043 147.6607838596043 161.5101241729302 173.8214905237701 184.6562641923786 + 234.5147959861216 254.6474401194407 266.4385386783982 275.1224493134379 282.1192618430768 + 144.3034887320350 147.0869026517344 147.8329014568262 147.9104560214784 147.6607838596043 +$TABLE_3 + 5 5 + 240.0000000000000 280.0000000000000 320.0000000000000 360.0000000000000 400.0000000000000 + 55000.00000000000 141250.0000000000 227500.0000000000 313750.0000000000 400000.0000000000 + 0.3481914037990427 0.4095738848363432 0.4702915495291079 0.5305832630019515 0.5905860753062941 + 0.1402248442644343 0.1561641066947830 0.1807449135220719 0.2048688905909286 0.2286916876895624 + 8.9273098371410883E-02 444320344302E-02 0.1107157321741076 0.1261131151219874 0.1411960889105597 + 6.5633087958005529E-02 6.7151770525718760E-02 7.9166190580894386E-02 9.0648541948697259E-02 0.1018021936800781 + 5.1900452073168755E-02 5.1900452073168755E-02 6.1203439903817498E-02 7.0470885355410795E-02 7.9394150509035255E-02 + 234.5147959861216 254.6474401194407 266.4385386783982 275.1224493134379 282.1192618430768 + 0.3397008120139503 0.1402248442644343 8.9273098371410883E-02 6.5633087958005529E-02 5.1900452073168755E-02 +$TABLE_4 + 5 5 + 240.0000000000000 280.0000000000000 320.0000000000000 + 360.0000000000000 400.0000000000000 + 55000.00000000000 141250.0000000000 227500.0000000000 + 313750.0000000000 400000.0000000000 + 644.8541874954598 719.7151797641214 790.5095338665236 + 857.2366695715677 919.8970271275263 674.7791960407488 + 721.6690087704523 792.3443057727616 858.9628531024553 + 921.5203908806476 698.8825060917976 723.6988581975571 + 794.2197800381258 860.7115954107751 923.1563593500967 + 716.9720727393429 725.8148416296926 796.1389707941393 + 862.4838781958006 924.8052629190319 731.8370306406068 + 731.8370306406068 798.1052760769886 864.2807537481988 + 926.4674455748727 + 234.5147959861216 254.6474401194407 266.4385386783982 + 275.1224493134379 282.1192618430768 + 634.2715599484291 674.7791960407488 698.8825060917976 + 716.9720727393429 731.8370306406068 +$TABLE_5 + 5 5 + 240.0000000000000 280.0000000000000 320.0000000000000 + 360.0000000000000 400.0000000000000 + 55000.00000000000 141250.0000000000 227500.0000000000 + 313750.0000000000 400000.0000000000 + 733.5432454120784 806.0266427729431 875.3572680256013 + 941.1190584653727 1003.111496904734 773.6370689985604 + 816.3711389524509 882.8108960718668 946.7561020285661 + 1007.512627818587 807.0628915610768 827.9980795292242 + 890.7840080089826 952.6232673331618 1012.020341520216 + 834.4867573220491 841.2008081815245 899.3408892587429 + 958.7368554939568 1016.639020749522 858.8621950799120 + 858.8621950799120 908.5575294975291 965.1148207606369 + 1021.373299533599 + 234.5147959861216 254.6474401194407 266.4385386783982 + 275.1224493134379 282.1192618430768 + 723.4007165554241 773.6370689985604 807.0628915610768 + 834.4867573220491 858.8621950799120 +$TABLE_6 + 5 5 + 240.0000000000000 280.0000000000000 320.0000000000000 + 360.0000000000000 400.0000000000000 + 55000.00000000000 141250.0000000000 227500.0000000000 + 313750.0000000000 400000.0000000000 + -154611.9428101669 -132557.4464307218 -115997.3804840207 + -103111.8528686964 -92801.69633261529 -959671.2544114459 + -873595.5996881851 -764854.3591592711 -680011.5179211511 + -612056.5815781737 -2374639.983863590 -2262401.182121443 + -1982996.713344134 -1763666.875581090 -1587621.548245155 + -4363498.860521937 -4291069.408039705 -3768294.493979145 + -3353439.511168068 -3019297.561697437 -6897297.389986640 + -6897297.389986640 -6117240.959939656 -5448329.203176694 + -4906781.742740840 + 234.5147959861216 254.6474401194407 266.4385386783982 + 275.1224493134379 282.1192618430768 + -158218.3298083634 -959671.2544114459 -2374639.983863590 + -4363498.860521937 -6897297.389986640 +$TABLE_7 + 5 5 + 240.0000000000000 280.0000000000000 320.0000000000000 + 360.0000000000000 400.0000000000000 + 55000.00000000000 141250.0000000000 227500.0000000000 + 313750.0000000000 400000.0000000000 + 1786.030158361605 1904.585719752846 2016.779512855124 + 2123.714780772325 2226.113661461630 1748.483241728080 + 1823.905112007991 1937.275071424368 2044.975937996867 + 2147.900948384664 1740.481415370240 1781.062492767768 + 1895.721351626920 2004.229786438356 2107.697911720815 + 1735.932309549977 1750.655644805153 1866.736275194941 + 1976.098388594144 2080.127578672003 1732.877658407725 + 1732.877658407725 1844.072019825660 1954.338672102516 + 2058.947921108171 + 234.5147959861216 254.6474401194407 266.4385386783982 + 275.1224493134379 282.1192618430768 + 1769.188083636742 1748.483241728080 1740.481415370240 + 1735.932309549977 1732.877658407725 +$TABLE_8 + 5 5 + 240.0000000000000 280.0000000000000 320.0000000000000 + 360.0000000000000 400.0000000000000 + 55000.00000000000 141250.0000000000 227500.0000000000 + 313750.0000000000 400000.0000000000 + 9.4859076076546339E-06 1.1098072605476093E-05 1.2682727526284134E-05 + 1.4225406873221409E-05 1.5718844908709814E-05 1.0113426130489583E-05 + 1.1128384712488209E-05 1.2707428862857411E-05 1.4246015436604838E-05 + 1.5736277818145713E-05 1.0622849104772714E-05 1.1161623164463395E-05 + 1.2733802630234965E-05 1.4267707066828921E-05 1.5754472822231917E-05 + 1.1006164112931968E-05 1.1198285735280963E-05 1.2761997824348757E-05 + 1.4290540551654799E-05 1.5773456592676736E-05 1.1321670531250205E-05 + 1.1321670531250205E-05 1.2792178337213511E-05 1.4314589371365448E-05 + 1.5793272119540815E-05 + 234.5147959861216 254.6474401194407 266.4385386783982 + 275.1224493134379 282.1192618430768 + 9.2638729380455867E-06 1.0113426130489583E-05 1.0622849104772714E-05 + 1.1006164112931968E-05 1.1321670531250205E-05 +$TABLE_9 + 5 5 + 240.0000000000000 280.0000000000000 320.0000000000000 + 360.0000000000000 400.0000000000000 + 55000.00000000000 141250.0000000000 227500.0000000000 + 313750.0000000000 400000.0000000000 + 8.7318132315962922E-03 1.2367010749671227E-02 1.5479028969457832E-02 + 1.8065003454933646E-02 2.0123880410625419E-02 1.0144646731520802E-02 + 1.2375365154798446E-02 1.5478277206731493E-02 1.8061024664563677E-02 + 2.0119085407073307E-02 1.1241534315205219E-02 1.2399427073087721E-02 + 1.5485543810413213E-02 1.8061449100116474E-02 2.0116675809444373E-02 + 1.2033532782386015E-02 1.2438788590707257E-02 1.5500573740125312E-02 + 1.8066241345977847E-02 2.0116840814255434E-02 1.2665788438301477E-02 + 1.2665788438301477E-02 1.5523214978202122E-02 1.8075268343013946E-02 + 2.0119597534766078E-02 + 234.5147959861216 254.6474401194407 266.4385386783982 + 275.1224493134379 282.1192618430768 + 8.1928902502712009E-03 1.0144646731520802E-02 1.1241534315205219E-02 + 1.2033532782386015E-02 1.2665788438301477E-02 +$$SAT_TABLE + 5 4 9 + 72429.02057212396 669962.0863334293 1267495.152094735 1865028.217856040 2462561.283617346 + 239.9999762675480 298.3735220799146 321.5184369582214 337.4761886205648 349.9999999385916 + 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 + 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 + 157792.0768617229 234423.4028017849 268281.8836128900 293391.1384832433 314607.4539245256 + 1226.028409346856 1411.621580596332 1538.623093445163 1688.584812454093 1901.162632743923 + 1394.714486693549 1204.132429873541 1105.717252974825 1023.909256983076 947.2617335591748 + 7.3133047353247582E-06 1.5554466347974513E-05 2.3546196865974526E-05 3.4054173679783401E-05 4.9653987218547186E-05 + 836.0313698668695 1119.425611702024 1226.996801517600 1301.491841955899 1361.442808503099 + 1087.583619605917 1180.740355522880 1205.127882782043 1215.472691163284 1218.190684502114 + 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 + 4.5254108276952397E-04 2.0762809729902491E-04 1.5172962190515123E-04 1.1920087721741867E-04 9.5575862436921458E-05 + 0.1106897637097442 8.1042418175005723E-02 6.9512380598733353E-02 6.1769627812551062E-02 5.5819411652293066E-02 + 380136.7395374610 415734.4668987678 426586.7644804174 431942.8510520858 434298.6702392859 + 736.4416185000921 926.5982109988553 1070.936300636404 1242.721436965689 1480.769140200031 + 3.808368423257648 32.02662176918217 61.85130414931603 94.90868505077478 132.4954540908795 + 5.4108323518481296E-05 \ No newline at end of file