diff --git a/test_problems/cathermo/HMW_graph_VvT/.cvsignore b/test_problems/cathermo/HMW_graph_VvT/.cvsignore new file mode 100644 index 000000000..da33a5c43 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_VvT/.cvsignore @@ -0,0 +1,11 @@ +Makefile +output.txt +outputa.txt +.cvsignore.swp +.depends +HMW_graph_VvT +HMW_graph_VvT.d +V_standalone +diff_test.out +sortAlgorithms.d +table.csv diff --git a/test_problems/cathermo/HMW_graph_VvT/HMW_NaCl_sp1977_alt.xml b/test_problems/cathermo/HMW_graph_VvT/HMW_NaCl_sp1977_alt.xml new file mode 100644 index 000000000..1ca6c066a --- /dev/null +++ b/test_problems/cathermo/HMW_graph_VvT/HMW_NaCl_sp1977_alt.xml @@ -0,0 +1,243 @@ + + + + + + H2O(L) Cl- H+ Na+ OH- + + + 298.15 + 101325.0 + + Na+:6.0954 + Cl-:6.0954 + H+:2.1628E-9 + OH-:1.3977E-6 + + + + + + + + + + + + 0.0765, 0.008946, -3.3158E-6, + -777.03, -4.4706 + + 0.2664, 6.1608E-5, 1.0715E-6 + 0.0 + 0.00127, -4.655E-5, 0.0, + 33.317, 0.09421 + + 2.0 + + + + 0.1775, 0.0, 0.0, 0.0, 0.0 + 0.2945, 0.0, 0.0 + 0.0 + 0.0008, 0.0, 0.0, 0.0, 0.0 + 2.0 + + + + 0.0864, 0.0, 0.0, 0.0, 0.0 + 0.253, 0.0, 0.0 + 0.0 + 0.0044, 0.0, 0.0, 0.0, 0.0 + 2.0 + + + + -0.05 + + + + -0.05 + -0.006 + + + + 0.036 + + + + 0.036 + -0.004 + + + + H2O(L) + + O H C Fe Si N Na Cl + + + + + + + + + + H:2 O:1 + + + + 7.255750050E+01, -6.624454020E-01, 2.561987460E-03, -4.365919230E-06, + 2.781789810E-09, -4.188654990E+04, -2.882801370E+02 + + + + + + 0.018068 + + + + + + Na:1 + +1 + + + + -57993.47558 , 305112.6040 , -592222.1591 , + 401977.9827 , 804.4195980 , 10625.24901 , + -133796.2298 + + + + + + + 0.00834 + + + + + + Cl:1 + -1 + + + + 0.00834 + + + + + 56696.2042 , -297835.978 , 581426.549 , + -401759.991 , -804.301136 , -10873.8257 , + 130650.697 + + + + + + + + H:1 + +1 + + 0.0 + + + + 0.0 + 3 + + 0.0 , 0.0, 0.0 + + + 273.15, 298.15 , 623.15 + + + + + + + + O:1 H:1 + -1 + + + 0.00834 + + + + + 44674.99961 , -234943.0414 , 460522.8260 , + -320695.1836 , -638.5044716 , -8683.955813 , + 102874.2667 + + + + + + + + diff --git a/test_problems/cathermo/HMW_graph_VvT/HMW_graph_VvT.cpp b/test_problems/cathermo/HMW_graph_VvT/HMW_graph_VvT.cpp new file mode 100755 index 000000000..fd1053aa6 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_VvT/HMW_graph_VvT.cpp @@ -0,0 +1,348 @@ +/** + * @file HMW_graph_VvT + */ + +/* + * $Author$ + * $Date$ + * $Revision$ + */ +#include + +#ifdef SRCDIRTREE +#include "ct_defs.h" +#include "logger.h" +#include "TemperatureTable.h" +#include "HMWSoln.h" +#include "importCTML.h" +#else +#include "vcsc.h" +#include "cantera/Cantera.h" +#include "vcs_Cantera_input.h" +#include "vcs_Cantera_convert.h" +#include "cantera/kernel/logger.h" +#include "cantera/thermo.h" +#include "TemperatureTable.h" +#include "ThermoPhase.h" +#include "HMWSoln.h" +#include "importCTML.h" + +#endif +using namespace Cantera; + +class fileLog: public Logger { +public: + fileLog(string fName) { + m_fName = fName; + m_fs.open(fName.c_str()); + } + + virtual void write(const string& msg) { + m_fs << msg; + m_fs.flush(); + } + + virtual ~fileLog() { + m_fs.close(); + } + + string m_fName; + ofstream m_fs; + +}; + +void printUsage() { + cout << "usage: HMW_test " << endl; + cout <<" -> Everything is hardwired" << endl; +} + +void pAtable(HMWSoln *HMW) { + int nsp = HMW->nSpecies(); + double acMol[100]; + double mf[100]; + double activities[100]; + double moll[100]; + + HMW->getMolalityActivityCoefficients(acMol); + HMW->getMoleFractions(mf); + HMW->getActivities(activities); + HMW->getMolalities(moll); + string sName; + printf(" Name Activity ActCoeffMolal " + " MoleFract Molality\n"); + for (int k = 0; k < nsp; k++) { + sName = HMW->speciesName(k); + printf("%16s %13g %13g %13g %13g\n", + sName.c_str(), activities[k], acMol[k], mf[k], moll[k]); + } +} + +int main(int argc, char **argv) +{ + + int retn = 0; + int i; + + try { + + char iFile[80]; + strcpy(iFile, "HMW_NaCl.xml"); + if (argc > 1) { + strcpy(iFile, argv[1]); + } + double V0[20], pmV[20]; + + //fileLog *fl = new fileLog("HMW_graph_1.log"); + //setLogger(fl); + + HMWSoln *HMW = new HMWSoln(iFile, "NaCl_electrolyte"); + + + /* + * Load in and initialize the + */ + Cantera::ThermoPhase *solid = newPhase("NaCl_Solid.xml","NaCl(S)"); + + + int nsp = HMW->nSpecies(); + //double acMol[100]; + //double act[100]; + double mf[100]; + double moll[100]; + HMW->getMoleFractions(mf); + string sName; + + TemperatureTable TTable(15, false, 273.15, 25., 0, 0); + + + HMW->setState_TP(298.15, 1.01325E5); + + int i1 = HMW->speciesIndex("Na+"); + int i2 = HMW->speciesIndex("Cl-"); + //int i3 = HMW->speciesIndex("H2O(L)"); + for (i = 0; i < nsp; i++) { + moll[i] = 0.0; + } + HMW->setMolalities(moll); + + + double ISQRT; + double Is = 0.0; + + /* + * Set the Pressure + */ + double pres = OneAtm; + + /* + * Fix the molality + */ + Is = 6.146; + ISQRT = sqrt(Is); + moll[i1] = Is; + moll[i2] = Is; + HMW->setState_TPM(298.15, pres, moll); + double Xmol[30]; + HMW->getMoleFractions(Xmol); + double meanMW = HMW->meanMolecularWeight(); + + /* + * ThermoUnknowns + */ + double T; + + double V0_NaCl, V0_Naplus, V0_Clminus, Delta_V0s, V0_H2O; + double V_NaCl, V_Naplus, V_Clminus, V_H2O; + double molarV0; +#ifdef DEBUG_HKM + FILE *ttt = fopen("table.csv","w"); +#endif + printf("A_V : Comparison to Pitzer's book, p. 99, can be made.\n"); + printf(" Agreement to 3 sig digits \n"); + printf("\n"); + + printf("Delta_V0: Heat Capacity of Solution per mole of salt (standard states)\n"); + printf(" rxn for the ss heat of soln: " + "NaCl(s) -> Na+(aq) + Cl-(aq)\n"); + + printf("\n"); + printf("Delta_Vs: Delta volume of Solution per mole of salt\n"); + printf(" rxn for heat of soln: " + " n1 H2O(l,pure) + n2 NaCl(s) -> n2 MX(aq) + n1 H2O(l) \n"); + printf(" Delta_Hs = (n1 h_H2O_bar + n2 h_MX_bar " + "- n1 h_H2O_0 - n2 h_MX_0)/n2\n"); + printf("\n"); + printf("phiV: phiV, calculated from the program, is checked\n"); + printf(" against analytical formula in V_standalone program.\n"); + printf(" (comparison against Pitzer book, p. 97, eqn. 96)\n"); + + /* + * Create a Table of NaCl Enthalpy Properties as a Function + * of the Temperature + */ + printf("\n\n"); + printf(" T, Pres, Aphi, A_V," + " Delta_V0," + " Delta_Vs, Vex, phiV," + " MolarV, MolarV0\n"); + printf(" Kelvin, bar, sqrt(kg/gmol),sqrt(kg/gmol)cm3/gmol," + "cm**3/gmolSalt," + "cm**3/gmolSalt,cm**3/gmolSoln,cm**3/gmolSalt," + "cm**3/gmol, cm**3/gmol\n"); +#ifdef DEBUG_HKM + fprintf(ttt,"T, Pres, A_V, Vex, phiV, MolarV, MolarV0\n"); + fprintf(ttt,"Kelvin, bar, sqrt(kg/gmol)cm3/gmol, cm3/gmolSoln, cm3/gmolSalt, kJ/gmolSoln," + "kJ/gmolSoln\n"); +#endif + for (i = 0; i < TTable.NPoints + 1; i++) { + if (i == TTable.NPoints) { + T = 323.15; + } else { + T = TTable.T[i]; + } + /* + * RT is in units of J/kmolK + */ + //double RT = GasConstant * T; + + /* + * Make sure we are at the saturation pressure or above. + */ + + double psat = HMW->satPressure(T); + + pres = OneAtm; + if (psat > pres) pres = psat; + + + HMW->setState_TPM(T, pres, moll); + + solid->setState_TP(T, pres); + + /* + * Get the Standard State volumes m3/kmol + */ + solid->getStandardVolumes(V0); + V0_NaCl = V0[0]; + HMW->getStandardVolumes(V0); + V0_H2O = V0[0]; + V0_Naplus = V0[i1]; + V0_Clminus = V0[i2]; + + /* + * Calculate the standard state volume change of solution + * for NaCl(s) -> Na+ + Cl- + * units: m3 / kmol + */ + Delta_V0s = V0_Naplus + V0_Clminus - V0_NaCl; + + double dd = solid->density(); + double MW_NaCl = solid->meanMolecularWeight(); + V_NaCl = MW_NaCl / dd; + //printf("V_NaCl = %g , V0_NaCl = %g %g\n", V_NaCl, V0_NaCl, 1.0/solid->molarDensity()); + + /* + * Get the partial molar volumes + */ + HMW->getPartialMolarVolumes(pmV); + V_H2O = pmV[0]; + V_Naplus = pmV[i1]; + V_Clminus = pmV[i2]; + + + //double Delta_V_Salt = V_NaCl - (V_Naplus + V_Clminus); + + /* + * Calculate the molar volume of solution + */ + double dsoln = HMW->density(); + meanMW = HMW->meanMolecularWeight(); + double molarV = meanMW / dsoln; + //double md = HMW->molarDensity(); + //printf("compare %g %g\n", molarV, 1.0/md); + + /* + * Calculate the delta volume of solution for the reaction + * NaCl(s) -> Na+ + Cl- + */ + double Delta_Vs = (Xmol[0] * V_H2O + + Xmol[i1] * V_Naplus + + Xmol[i2] * V_Clminus + - Xmol[0] * V0_H2O + - Xmol[i1] * V_NaCl); + Delta_Vs /= Xmol[i1]; + + + /* + * Calculate the apparent molar volume, J, from the + * partial molar quantities, units m3/kmol + */ + double Vex = (Xmol[0] * (V_H2O - V0_H2O) + + Xmol[i1] * (V_Naplus - V0_Naplus) + + Xmol[i2] * (V_Clminus - V0_Clminus)); + + /* + * Calculate the apparent relative molal volume, phiV, + * units of m3/kmol + */ + double phiV = Vex / Xmol[i1]; + + double Aphi = HMW->A_Debye_TP(T, pres) / 3.0; + //double AL = HMW->ADebye_L(T,pres); + double Av = HMW->ADebye_V(T, pres) * 1.0E3; + + + molarV0 = 0.0; + for (int k = 0; k < nsp; k++) { + molarV0 += Xmol[k] * V0[k]; + } + + if (i != TTable.NPoints+1) { + printf("%13g, %13g, %13g, %13g, %13g, %13g, " + "%13g, %13g, %13g, %13g\n", + T, pres*1.0E-5, Aphi, Av, Delta_V0s*1.0E3, Delta_Vs*1.0E3, + Vex*1.0E3, phiV*1.0E3, molarV*1.0E3 , molarV0*1.0E3 ); +#ifdef DEBUG_HKM + fprintf(ttt,"%g, %g, %g, %g, %g, %g, %g\n", + T, pres*1.0E-5, Av, Vex*1.0E3, phiV*1.0E3, molarV*1.0E3 , molarV0*1.0E3); +#endif + } + + } + + printf("Breakdown of Volume Calculation at 323.15 K, 1atm:\n"); + + printf(" Species MoleFrac Molal V0 " + " partV (partV - V0)\n"); + printf(" H2O(L)"); + printf("%13g %13g %13g %13g %13g\n", Xmol[0], moll[0], V0_H2O*1.E3, V_H2O*1.E3, + (V_H2O-V0_H2O)*1.E3); + printf(" Na+ "); + printf("%13g %13g %13g %13g %13g\n", Xmol[i1], moll[i1], + V0_Naplus*1.E3 , V_Naplus*1.E3, (V_Naplus -V0_Naplus)*1.E3); + printf(" Cl- "); + printf("%13g %13g %13g %13g %13g\n", Xmol[i2], moll[i2], + V0_Clminus*1.E3, V_Clminus*1.E3, (V_Clminus - V0_Clminus)*1.E3); + + printf(" NaCl(s)"); + printf("%13g %13g %13g %13g\n", 1.0, + V0_NaCl*1.E3 , V_NaCl*1.E3, V_NaCl*1.E3 - V0_NaCl*1.E3); + + + delete HMW; + HMW = 0; + delete solid; + solid = 0; + Cantera::appdelete(); +#ifdef DEBUG_HKM + fclose(ttt); +#endif + return retn; + + } catch (CanteraError) { + printf("caught error\n"); + showErrors(); + Cantera::appdelete(); + return -1; + } +} diff --git a/test_problems/cathermo/HMW_graph_VvT/Makefile.in b/test_problems/cathermo/HMW_graph_VvT/Makefile.in new file mode 100644 index 000000000..506765552 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_VvT/Makefile.in @@ -0,0 +1,116 @@ +#!/bin/sh + +############################################################################ +# +# Makefile to compile and link a C++ application to +# Cantera. +# +############################################################################# + +# addition to suffixes +.SUFFIXES : .d + +# the name of the executable program to be created +PROG_NAME = HMW_graph_VvT + +# the object files to be linked together. List those generated from Fortran +# and from C/C++ separately +OBJS = HMW_graph_VvT.o sortAlgorithms.o + +# Location of the current build. Will assume that tests are run +# in the source directory tree location +src_dir_tree = 1 + +# additional flags to be passed to the linker. If your program +# requires other external libraries, put them here +LINK_OPTIONS = @EXTRA_LINK@ + +############################################################################# + +# Check to see whether we are in the msvc++ environment +os_is_win = @OS_IS_WIN@ + +# Fortran libraries +FORT_LIBS = @FLIBS@ + +# the C++ compiler +CXX = @CXX@ + +# C++ compile flags +ifeq ($(src_dir_tree), 1) +CXX_FLAGS = -DSRCDIRTREE @CXXFLAGS@ +else +CXX_FLAGS = @CXXFLAGS@ +endif + +# Ending C++ linking libraries +LCXX_END_LIBS = @LCXX_END_LIBS@ + +# the directory where the Cantera libraries are located +CANTERA_LIBDIR=@buildlib@ + +# required Cantera libraries +CANTERA_LIBS = @LOCAL_LIBS@ -lctcxx + +# the directory where Cantera include files may be found. +ifeq ($(src_dir_tree), 1) +CANTERA_INCDIR=../../../Cantera/src +INCLUDES=-I$(CANTERA_INCDIR) -I$(CANTERA_INCDIR)/thermo +else +CANTERA_INCDIR=@ctroot@/build/include/cantera +INCLUDES=-I$(CANTERA_INCDIR) -I$(CANTERA_INCDIR)/kernel +endif + +# flags passed to the C++ compiler/linker for the linking step +LCXX_FLAGS = -L$(CANTERA_LIBDIR) @LOCAL_LIB_DIRS@ @CXXFLAGS@ + +# How to compile C++ source files to object files +.@CXX_EXT@.@OBJ_EXT@: + $(CXX) -c $< $(INCLUDES) $(CXX_FLAGS) + +# How to compile the dependency file +.cpp.d: + g++ -MM $(INCLUDES) $(CXX_FLAGS) $*.cpp > $*.d + +# List of dependency files to be created +DEPENDS=$(OBJS:.o=.d) + +# Program Name +PROGRAM = $(PROG_NAME)$(EXE_EXT) + +all: $(PROGRAM) .depends V_standalone + +$(PROGRAM): $(OBJS) $(CANTERA_LIBDIR)/libcantera.a \ + $(CANTERA_LIBDIR)/libcaThermo.a + $(CXX) -o $(PROGRAM) $(OBJS) $(LCXX_FLAGS) $(LINK_OPTIONS) \ + $(CANTERA_LIBS) @LIBS@ $(FORT_LIBS) \ + $(LCXX_END_LIBS) + +V_standalone: V_standalone.o + $(CXX) -o V_standalone V_standalone.o \ + $(LCXX_FLAGS) $(LINK_OPTIONS) $(LCXX_END_LIBS) + +# depends target -> forces recalculation of dependencies +depends: + @MAKE@ .depends + +.depends: $(DEPENDS) + cat $(DEPENDS) > .depends + +# Do the test -> For the windows vc++ environment, we have to skip checking on +# whether the program is uptodate, because we don't utilize make +# in that environment to build programs. +test: +ifeq ($(os_is_win), 1) +else + @MAKE@ $(PROGRAM) +endif + ./runtest + +clean: + $(RM) $(OBJS) $(PROGRAM) $(DEPENDS) .depends + ../../../bin/rm_cvsignore + (if test -d SunWS_cache ; then \ + $(RM) -rf SunWS_cache ; \ + fi ) + diff --git a/test_problems/cathermo/HMW_graph_VvT/NaCl_Solid.xml b/test_problems/cathermo/HMW_graph_VvT/NaCl_Solid.xml new file mode 100644 index 000000000..d711be8ff --- /dev/null +++ b/test_problems/cathermo/HMW_graph_VvT/NaCl_Solid.xml @@ -0,0 +1,39 @@ + + + + + + + + + O H C Fe Ca N Na Cl + + NaCl(S) + + 2.165 + + + + + + + + + + + Na:1 Cl:1 + + + + 50.72389, 6.672267, -2.517167, + 10.15934, -0.200675, -427.2115, + 130.3973 + + + + 2.165 + + + + + diff --git a/test_problems/cathermo/HMW_graph_VvT/README b/test_problems/cathermo/HMW_graph_VvT/README new file mode 100644 index 000000000..f22bf8f5e --- /dev/null +++ b/test_problems/cathermo/HMW_graph_VvT/README @@ -0,0 +1 @@ +Check on the enthalpy routines for HMWSoln diff --git a/test_problems/cathermo/HMW_graph_VvT/TemperatureTable.h b/test_problems/cathermo/HMW_graph_VvT/TemperatureTable.h new file mode 100644 index 000000000..70cb2bf56 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_VvT/TemperatureTable.h @@ -0,0 +1,129 @@ +/* + * $Id$ + */ +/* + * Copywrite 2004 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000, there is a non-exclusive license for use of this + * work by or on behalf of the U.S. Government. Export of this program + * may require a license from the United States Government. + */ + +#ifndef TEMPERATURE_TABLE_H +#define TEMPERATURE_TABLE_H +#include "sortAlgorithms.h" +//#include "mdp_allo.h" +#include +using std::vector; + +/***********************************************************************/ +/***********************************************************************/ +/***********************************************************************/ +/** + * This Class constructs a vector of temperature from which to make + * a table. + */ +class TemperatureTable { + +public: + int NPoints; + bool Include298; + double Tlow; //!< Min temperature for thermo data fit + double Thigh; //!< Max temperature for thermo table + double DeltaT; + vector T; + int numAddedTs; + vector AddedTempVector; +public: + /* + * Default constructor for TemperatureTable() + */ + TemperatureTable(const int nPts = 14, + const bool inc298 = true, + const double tlow = 300., + const double deltaT = 100., + const int numAdded = 0, + const double *addedTempVector = 0) : + NPoints(nPts), + Include298(inc298), + Tlow(tlow), + DeltaT(deltaT), + T(0), + numAddedTs(numAdded) { + /****************************/ + int i; + // AddedTempVector = mdp_alloc_dbl_1(numAdded, 0.0); + AddedTempVector.resize(numAdded, 0.0); + for (int i = 0; i < numAdded; i++) { + AddedTempVector[i] = addedTempVector[i]; + } + //mdp_copy_dbl_1(AddedTempVector, addedTempVector, numAdded); + // T = mdp_alloc_dbl_1(NPoints, 0.0); + T.resize(NPoints, 0.0); + double TCurrent = Tlow; + for (i = 0; i < NPoints; i++) { + T[i] = TCurrent; + TCurrent += DeltaT; + } + if (Include298) { + T.push_back(298.15); + //mdp_realloc_dbl_1(&T, NPoints+1, NPoints, 298.15); + NPoints++; + } + if (numAdded > 0) { + //mdp_realloc_dbl_1(&T, NPoints+numAdded, NPoints, 0.0); + T.resize( NPoints+numAdded, 0.0); + for (i = 0; i < numAdded; i++) { + T[i+NPoints] = addedTempVector[i]; + } + NPoints += numAdded; + } + + sort_dbl_1(DATA_PTR(T), NPoints); + + + } + /***********************************************************************/ + /***********************************************************************/ + /***********************************************************************/ + /* + * Destructor + */ + ~TemperatureTable() { + //mdp_safe_free((void **) &AddedTempVector); + // mdp_safe_free((void **) &T); + } + + /***********************************************************************/ + /***********************************************************************/ + /***********************************************************************/ + /* + * Overloaded operator[] + * + * return the array value in the vector + */ + double operator[](const int i) { + return T[i]; + } + /***********************************************************************/ + /***********************************************************************/ + /***********************************************************************/ + /* + * size() + */ + int size() { + return NPoints; + } +/***********************************************************************/ +/***********************************************************************/ +/***********************************************************************/ + /* + * Block assignment and copy constructors: not needed. + */ +private: + TemperatureTable(const TemperatureTable &); + TemperatureTable& operator=(const TemperatureTable&); +}; +/***********************************************************************/ +/***********************************************************************/ +/***********************************************************************/ +#endif diff --git a/test_problems/cathermo/HMW_graph_VvT/V_standalone.cpp b/test_problems/cathermo/HMW_graph_VvT/V_standalone.cpp new file mode 100644 index 000000000..bc319eb57 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_VvT/V_standalone.cpp @@ -0,0 +1,182 @@ + +#include +#include +#include + +using namespace std; + + + +/* + * Values of A_V : tabular form + * units sqrt(kg/gmol)cm3/gmol + */ +double A_V(double temp) { + double retn; + if (temp == 323.15) { + retn = 2.37356; + } else if (temp == 473.15) { + retn = 15.0766; + } else { + printf("A_V unknown temp value %g\n", temp); + exit(-1); + } + return retn; +} + +double Beta0(double temp, int ifunc) { + double q1 = 0.0765; + double q2 = -777.03; + double q3 = -4.4706; + double q4 = 0.008946; + double q5 = -3.3158E-6; + double retn; + double tref = 298.15; + if (ifunc == 0) { + retn = q1 + q2 * (1.0/temp - 1.0/tref) + + q3 * (log(temp/tref)) + q4 * (temp - tref) + + q5 * (temp * temp - tref * tref); + } else if (ifunc == 1) { + retn = (- q2 * 1.0/(temp* temp) + + q3 / temp + + q4 + + 2.0 * temp * q5); + } else if (ifunc == 2) { + retn = ( 2.0 * q2 * 1.0/(temp* temp*temp) + - q3 / (temp*temp) + + 2.0 * q5); + } else if (ifunc == 3) { + retn = 0.0; + } else { + exit(-1); + } + return retn; +} + +double Beta1(double temp, int ifunc) { + double q6 = 0.2664; + double q9 = 6.1608E-5; + double q10 = 1.0715E-6; + double retn; + double tref = 298.15; + if (ifunc == 0) { + retn = q6 + q9 * (temp - tref) + + q10 * (temp * temp - tref * tref); + } else if (ifunc == 1) { + retn = q9 + 2.0 * q10 * temp; + } else if (ifunc == 2) { + retn = 2.0 * q10; + } else if (ifunc == 3) { + retn = 0.0; + } else { + exit(-1); + } + return retn; +} + +double Cphi(double temp, int ifunc) { + double q11 = 0.00127; + double q12 = 33.317; + double q13 = 0.09421; + double q14 = -4.655E-5; + double retn; + double tref = 298.15; + if (ifunc == 0) { + retn = q11 + q12 * (1.0/temp - 1.0/tref) + + q13 * (log(temp/tref)) + q14 * (temp - tref); + } else if (ifunc == 1) { + retn = - q12 / (temp * temp) + + q13 / temp + q14; + } else if (ifunc == 2) { + retn = + 2.0 * q12 / (temp * temp * temp) + - q13 / (temp * temp) ; + } else if (ifunc == 3) { + retn = 0.0; + } else { + exit(-1); + } + return retn; +} + +double calc(double temp, double Iionic) { + /* + * Gas Constant in J gmol-1 K-1 + */ + double GasConst = 8.314472; + + double Aphi = 0.0; + if (temp == 323.15) { + Aphi = 0.4102995331359; + } else if (temp == 473.15) { + Aphi = 0.622777; + } else { + printf("ERROR: unknown temp\n"); + exit(-1); + } + + + /* + * Calculate A_V in sqrt(kg/gmol)cm3/gmol + */ + double Av = A_V(temp); + + double beta0prime3 = Beta0(temp, 3); + printf(" beta0prime3 = %g\n", beta0prime3); + + double beta1prime3 = Beta1(temp, 3); + printf(" beta1prime3 = %g\n", beta1prime3); + + double cphiprime3= Cphi(temp, 3); + printf(" Cphiprime = %g\n", cphiprime3); + + double vm = 1.0; + double vx = 1.0; + double v = vm + vx; + double m = Iionic; + double zm = 1.; + double zx = 1.0; + + double sqrtI = sqrt(Iionic); + + double alpha = 2.0; + double a2 = alpha * alpha; + double b = 1.2; + + double BVmx = beta0prime3 + 2.0 * beta1prime3 / (a2* Iionic) * + (1.0 - (1.0 + alpha * sqrtI) * exp(-alpha*sqrtI) ); + + double CVmx = 0.5 * sqrt(vm * vx) * cphiprime3; + + double phiV = v * zm * zx * (Av/(2.*b)) * log(1 + 1.2 * sqrtI) + + 2 * vm * vx * GasConst * temp *1.0E3 * ( m * BVmx + m * m * CVmx); + + printf(" phiV = %15.8g cm3/gmolSalt\n", phiV); + + double molecWeight = 18.01528; + + double RT = GasConst * temp * 1.0E-3; + + + double xo = 1.0 / (molecWeight/1000. * 2 * m + 1.0); + printf(" no = %g\n", xo); + + return phiV; +} + +main() { + + printf("Standalone test of the apparent relative molal excess volume, phiV:\n"); + printf(" (Check against simple formula in \n"); + printf(" Activity Coefficients in Eletrolyte Solutions, 2nd Ed K. S. Pitzer, " + "CRC Press, Boca Raton, 1991 \n"); + printf("T = 50C\n"); + double Iionic = 6.146; + printf("Ionic Strength = %g\n", Iionic); + + double res = calc(273.15 + 50., Iionic); + printf("T = 200C\n"); + printf("Ionic Strength = %g\n", Iionic); + + res = calc(273.15 + 200., Iionic); + +} diff --git a/test_problems/cathermo/HMW_graph_VvT/output_blessed.txt b/test_problems/cathermo/HMW_graph_VvT/output_blessed.txt new file mode 100644 index 000000000..b72613c52 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_VvT/output_blessed.txt @@ -0,0 +1,39 @@ +A_V : Comparison to Pitzer's book, p. 99, can be made. + Agreement to 3 sig digits + +Delta_V0: Heat Capacity of Solution per mole of salt (standard states) + rxn for the ss heat of soln: NaCl(s) -> Na+(aq) + Cl-(aq) + +Delta_Vs: Delta volume of Solution per mole of salt + rxn for heat of soln: n1 H2O(l,pure) + n2 NaCl(s) -> n2 MX(aq) + n1 H2O(l) + Delta_Hs = (n1 h_H2O_bar + n2 h_MX_bar - n1 h_H2O_0 - n2 h_MX_0)/n2 + +phiV: phiV, calculated from the program, is checked + against analytical formula in V_standalone program. + (comparison against Pitzer book, p. 97, eqn. 96) + + + T, Pres, Aphi, A_V, Delta_V0, Delta_Vs, Vex, phiV, MolarV, MolarV0 + Kelvin, bar, sqrt(kg/gmol),sqrt(kg/gmol)cm3/gmol,cm**3/gmolSalt,cm**3/gmolSalt,cm**3/gmolSoln,cm**3/gmolSalt,cm**3/gmol, cm**3/gmol + 273.15, 1.01325, 0.376717, 1.5062, -10.3142, -8.58207, 0.157016, 1.73214, 16.4205, 16.2635 + 298.15, 1.01325, 0.391447, 1.87435, -10.3142, -8.1587, 0.195394, 2.15551, 16.5003, 16.3049 + 323.15, 1.01325, 0.410293, 2.37352, -10.3142, -7.58465, 0.247431, 2.72956, 16.6872, 16.4398 + 348.15, 1.01325, 0.433273, 3.06983, -10.3142, -6.78389, 0.320018, 3.53032, 16.9618, 16.6418 + 373.15, 1.01418, 0.460559, 4.05167, -10.3142, -5.65477, 0.422371, 4.65945, 17.3246, 16.9022 + 398.15, 2.32238, 0.492454, 5.4513, -10.3142, -4.04518, 0.568278, 6.26903, 17.7872, 17.2189 + 423.15, 4.76165, 0.529514, 7.48096, -10.3142, -1.71106, 0.779862, 8.60316, 18.3759, 17.596 + 448.15, 8.92602, 0.572549, 10.4889, -10.3142, 1.74806, 1.09343, 12.0623, 19.1351, 18.0417 + 473.15, 15.5493, 0.622769, 15.0764, -10.3142, 7.02378, 1.57166, 17.338, 20.1415, 18.5698 + 498.15, 25.4972, 0.682036, 22.3389, -10.3142, 15.3757, 2.32875, 25.6899, 21.531, 19.2022 + 523.15, 39.7617, 0.753389, 34.4316, -10.3142, 29.2824, 3.58936, 39.5966, 23.5634, 19.974 + 548.15, 59.4639, 0.842213, 56.0339, -10.3142, 54.1252, 5.84132, 64.4394, 26.7856, 20.9443 + 573.15, 85.879, 0.959258, 98.8138, -10.3142, 103.322, 10.301, 113.637, 32.5242, 22.2232 + 598.15, 120.51, 1.13023, 198.467, -10.3142, 217.925, 20.6895, 228.239, 44.7424, 24.053 + 623.15, 165.294, 1.43872, 494.888, -10.3142, 558.811, 51.5903, 569.125, 78.7661, 27.1758 + 323.15, 1.01325, 0.410293, 2.37352, -10.3142, -7.58465, 0.247431, 2.72956, 16.6872, 16.4398 +Breakdown of Volume Calculation at 323.15 K, 1atm: + Species MoleFrac Molal V0 partV (partV - V0) + H2O(L) 0.818703 0 18.2334 18.1515 -0.0819526 + Na+ 0.0906484 6.146 8.34 10.0749 1.73486 + Cl- 0.0906484 6.146 8.34 10.0749 1.73486 + NaCl(s) 1 26.9942 26.9942 0 diff --git a/test_problems/cathermo/HMW_graph_VvT/runtest b/test_problems/cathermo/HMW_graph_VvT/runtest new file mode 100755 index 000000000..7d2c4b3a2 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_VvT/runtest @@ -0,0 +1,42 @@ +#!/bin/sh +# +# +temp_success="1" +/bin/rm -f output.txt outputa.txt + +########################################################################## +prog=HMW_graph_VvT +if test ! -x $prog ; then + echo $prog ' does not exist' + exit -1 +fi +########################################################################## +/bin/rm -f test.out test.diff output.txt + +################################################################# +# +CANTERA_DATA=${CANTERA_DATA:=../../../data/inputs}; export CANTERA_DATA +CANTERA_BIN=${CANTERA_BIN:=../../../bin} + +################################################################# + +$prog HMW_NaCl_sp1977_alt.xml > output.txt +retnStat=$? +if [ $retnStat != "0" ] +then + temp_success="0" + echo "$prog returned with bad status, $retnStat, check output" +fi + +$CANTERA_BIN/exp3to2.sh output.txt > outputa.txt +diff -w outputa.txt output_blessed.txt > diff_test.out +retnStat=$? +if [ $retnStat = "0" ] +then + echo "successful diff comparison on $prog test" +else + echo "unsuccessful diff comparison on $prog test" + echo "FAILED" > csvCode.txt + temp_success="0" +fi + diff --git a/test_problems/cathermo/HMW_graph_VvT/sortAlgorithms.cpp b/test_problems/cathermo/HMW_graph_VvT/sortAlgorithms.cpp new file mode 100644 index 000000000..d97e51b40 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_VvT/sortAlgorithms.cpp @@ -0,0 +1,54 @@ +/* + * @file sortAlgorithms.h + * + * $Author$ + * $Revision$ + * $Date$ + */ +/* + * Copywrite 2004 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government + * retains certain rights in this software. + * See file License.txt for licensing information. + */ + +#include "sortAlgorithms.h" + +/**************************************************************/ + +void sort_dbl_1(double * const x, const int n) { + double rra; + int ll = n/2; + int iret = n - 1; + while (1 > 0) { + if (ll > 0) { + ll--; + rra = x[ll]; + } else { + rra = x[iret]; + x[iret] = x[0]; + iret--; + if (iret == 0) { + x[0] = rra; + return; + } + } + int i = ll; + int j = ll + ll + 1; + while (j <= iret) { + if (j < iret) { + if (x[j] < x[j+1]) + j++; + } + if (rra < x[j]) { + x[i] = x[j]; + i = j; + j = j + j + 1; + } else { + j = iret + 1; + } + } + x[i] = rra; + } +} +/*****************************************************/ diff --git a/test_problems/cathermo/HMW_graph_VvT/sortAlgorithms.h b/test_problems/cathermo/HMW_graph_VvT/sortAlgorithms.h new file mode 100644 index 000000000..72a7fc2a2 --- /dev/null +++ b/test_problems/cathermo/HMW_graph_VvT/sortAlgorithms.h @@ -0,0 +1,21 @@ +/* + * @file sortAlgorithms.h + * + * $Author$ + * $Revision$ + * $Date$ + */ +/* + * Copywrite 2004 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government + * retains certain rights in this software. + * See file License.txt for licensing information. + */ + +#ifndef SORTALGORITHMS_H +#define SORTALGORITHMS_H + + +void sort_dbl_1(double * const x, const int n); + +#endif