GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2025-02-25 18:55:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

Added the excess volume test.

Browse Source
This commit is contained in:
Harry Moffat
2006-07-06 21:12:53 +00:00
parent ce87b6434c
commit 480829a8a7
12 changed files with 1225 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
test_problems/cathermo/HMW_graph_VvT/.cvsignore Normal file
View File

@@ -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

243
test_problems/cathermo/HMW_graph_VvT/HMW_NaCl_sp1977_alt.xml Normal file
View File

@@ -0,0 +1,243 @@
<?xml version="1.0"?>
<!--
$Id$
$Source$
$Name$
NaCl modeling Based on the Silvester&Pitzer 1977 treatment:
(L. F. Silvester, K. S. Pitzer, "Thermodynamics of Electrolytes:
8. High-Temperature Properties, including Enthalpy and Heat
Capacity, with application to sodium chloride",
J. Phys. Chem., 81, 19 1822 - 1828 (1977)
This modification reworks the Na+ standard state shomate
polynomial, so that the resulting DeltaG0 for the NaCl(s) -> Na+ + Cl-
reaction agrees closely with Silvester and Pitzer. The main
effect that this has is to change the predicted Na+ heat capacity
at low temperatures.
-->
<ctml>
<phase id="NaCl_electrolyte" dim="3">
<speciesArray datasrc="#species_waterSolution">
H2O(L) Cl- H+ Na+ OH-
</speciesArray>
<state>
<temperature units="K"> 298.15 </temperature>
<pressure units="Pa"> 101325.0 </pressure>
<soluteMolalities>
Na+:6.0954
Cl-:6.0954
H+:2.1628E-9
OH-:1.3977E-6
</soluteMolalities>
</state>
<thermo model="HMW">
<standardConc model="solvent_volume" />
<activityCoefficients model="Pitzer" TempModel="complex1">
<!-- Pitzer Coefficients
These coefficients are from Pitzer's main
paper, in his book.
-->
<A_Debye model="water" />
<ionicRadius default="3.042843" units="Angstroms">
</ionicRadius>
<binarySaltParameters cation="Na+" anion="Cl-">
<beta0> 0.0765, 0.008946, -3.3158E-6,
-777.03, -4.4706
</beta0>
<beta1> 0.2664, 6.1608E-5, 1.0715E-6 </beta1>
<beta2> 0.0 </beta2>
<Cphi> 0.00127, -4.655E-5, 0.0,
33.317, 0.09421
</Cphi>
<Alpha1> 2.0 </Alpha1>
</binarySaltParameters>
<binarySaltParameters cation="H+" anion="Cl-">
<beta0> 0.1775, 0.0, 0.0, 0.0, 0.0</beta0>
<beta1> 0.2945, 0.0, 0.0 </beta1>
<beta2> 0.0 </beta2>
<Cphi> 0.0008, 0.0, 0.0, 0.0, 0.0 </Cphi>
<Alpha1> 2.0 </Alpha1>
</binarySaltParameters>
<binarySaltParameters cation="Na+" anion="OH-">
<beta0> 0.0864, 0.0, 0.0, 0.0, 0.0 </beta0>
<beta1> 0.253, 0.0, 0.0 </beta1>
<beta2> 0.0 </beta2>
<Cphi> 0.0044, 0.0, 0.0, 0.0, 0.0 </Cphi>
<Alpha1> 2.0 </Alpha1>
</binarySaltParameters>
<thetaAnion anion1="Cl-" anion2="OH-">
<Theta> -0.05 </Theta>
</thetaAnion>
<psiCommonCation cation="Na+" anion1="Cl-" anion2="OH-">
<Theta> -0.05 </Theta>
<Psi> -0.006 </Psi>
</psiCommonCation>
<thetaCation cation1="Na+" cation2="H+">
<Theta> 0.036 </Theta>
</thetaCation>
<psiCommonAnion anion="Cl-" cation1="Na+" cation2="H+">
<Theta> 0.036 </Theta>
<Psi> -0.004 </Psi>
</psiCommonAnion>
</activityCoefficients>
<solvent> H2O(L) </solvent>
</thermo>
<elementArray datasrc="elements.xml"> O H C Fe Si N Na Cl </elementArray>
<kinetics model="none" >
</kinetics>
</phase>
<speciesData id="species_waterSolution">
<species name="H2O(L)">
<!-- H2O(L) liquid standard state -> pure H2O
The origin of the NASA polynomial is a bit murky. It does
fit the vapor pressure curve at 298K adequately.
-->
<atomArray>H:2 O:1 </atomArray>
<thermo>
<NASA Tmax="600.0" Tmin="273.14999999999998" P0="100000.0">
<floatArray name="coeffs" size="7">
7.255750050E+01, -6.624454020E-01, 2.561987460E-03, -4.365919230E-06,
2.781789810E-09, -4.188654990E+04, -2.882801370E+02
</floatArray>
</NASA>
</thermo>
<standardState model="waterPDSS">
<!--
Molar volume in m3 kmol-1.
(this is from Pitzer, Peiper, and Busey. However,
the result can be easily derived from ~ 1gm/cm**3)
-->
<molarVolume> 0.018068 </molarVolume>
</standardState>
</species>
<species name="Na+">
<!-- Na+ rework. Differences in the delta_G0 reaction
for salt formation were dumped into this polynomial.
-->
<atomArray> Na:1 </atomArray>
<charge> +1 </charge>
<thermo>
<Shomate Pref="1 bar" Tmax=" 593.15" Tmin=" 293.15">
<floatArray size="7">
-57993.47558 , 305112.6040 , -592222.1591 ,
401977.9827 , 804.4195980 , 10625.24901 ,
-133796.2298
</floatArray>
</Shomate>
</thermo>
<standardState model="constant_incompressible">
<!-- Na+ (aq) molar volume
Molar volume in m3 kmol-1.
(this is from Pitzer, Peiper, and Busey. We divide
NaCl (aq) value by 2 to get this)
-->
<molarVolume> 0.00834 </molarVolume>
</standardState>
</species>
<species name="Cl-">
<!-- Cl- (aq) standard state based on the unity molality convention
The shomate polynomial was created from the SUPCRT92
J. Phys Chem Ref article, and the codata recommended
values. DelHf(298.15) = -167.08 kJ/gmol
S(298.15) = 56.60 J/gmolK
There was a slight discrepency between those two, which was
resolved in favor of codata.
Notes: the order of the polynomials can be decreased by
dropping terms from the complete Shomate poly.
-->
<atomArray> Cl:1 </atomArray>
<charge> -1 </charge>
<standardState model="constant_incompressible">
<!-- Cl- (aq) molar volume
Molar volume in m3 kmol-1.
(this is from Pitzer, Peiper, and Busey. We divide
NaCl (aq) value by 2 to get this)
-->
<molarVolume> 0.00834 </molarVolume>
</standardState>
<thermo>
<Shomate Pref="1 atm" Tmax=" 623.15" Tmin=" 298.00">
<floatArray size="7">
56696.2042 , -297835.978 , 581426.549 ,
-401759.991 , -804.301136 , -10873.8257 ,
130650.697
</floatArray>
</Shomate>
</thermo>
</species>
<species name="H+">
<!-- H+ (aq) standard state based on the unity molality convention
The H+ standard state is set to zeroes by convention. This
includes it's contribution to the molar volume of solution.
-->
<atomArray> H:1 </atomArray>
<charge> +1 </charge>
<standardState model="constant_incompressible">
<molarVolume> 0.0 </molarVolume>
</standardState>
<thermo>
<Mu0 Pref="100000.0" Tmax="625.15." Tmin="273.15">
<H298 units="cal/mol"> 0.0 </H298>
<numPoints> 3 </numPoints>
<floatArray size="3" title="Mu0Values" units="Dimensionless">
0.0 , 0.0, 0.0
</floatArray>
<floatArray size="3" title="Mu0Temperatures">
273.15, 298.15 , 623.15
</floatArray>
</Mu0>
</thermo>
</species>
<species name="OH-">
<!-- OH- (aq) standard state based on the unity molality convention
The shomate polynomial was created with data from the SUPCRT92
J. Phys Chem Ref article, and from the codata recommended
values. DelHf(298.15) = -230.015 kJ/gmol
S(298.15) = -10.90 J/gmolK
There was a slight discrepency between those two, which was
resolved in favor of codata.
Notes: the order of the polynomials can be decreased by
dropping terms from the complete Shomate poly.
-->
<atomArray> O:1 H:1 </atomArray>
<charge> -1 </charge>
<standardState model="constant_incompressible">
<!-- OH- (aq) molar volume
This value is currently made up.
-->
<molarVolume> 0.00834 </molarVolume>
</standardState>
<thermo>
<Shomate Pref="1 atm" Tmax=" 623.15" Tmin=" 298.00">
<floatArray size="7">
44674.99961 , -234943.0414 , 460522.8260 ,
-320695.1836 , -638.5044716 , -8683.955813 ,
102874.2667
</floatArray>
</Shomate>
</thermo>
</species>
</speciesData>
</ctml>

348
test_problems/cathermo/HMW_graph_VvT/HMW_graph_VvT.cpp Executable file
View File

@@ -0,0 +1,348 @@
/**
* @file HMW_graph_VvT
*/
/*
* $Author$
* $Date$
* $Revision$
*/
#include <stdio.h>
#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;
}
}

116
test_problems/cathermo/HMW_graph_VvT/Makefile.in Normal file
View File

@@ -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 )

39
test_problems/cathermo/HMW_graph_VvT/NaCl_Solid.xml Normal file
View File

@@ -0,0 +1,39 @@
<?xml version="1.0"?>
<ctml>
<validate reactions="yes" species="yes"/>
<!-- phase NaCl(S) -->
<phase dim="3" id="NaCl(S)">
<elementArray datasrc="elements.xml">
O H C Fe Ca N Na Cl
</elementArray>
<speciesArray datasrc="#species_NaCl(S)"> NaCl(S) </speciesArray>
<thermo model="StoichSubstance">
<density units="g/cm3">2.165</density>
</thermo>
<transport model="None"/>
<kinetics model="none"/>
</phase>
<!-- species definitions -->
<speciesData id="species_NaCl(S)">
<!-- species NaCl(S) -->
<species name="NaCl(S)">
<atomArray> Na:1 Cl:1 </atomArray>
<thermo>
<Shomate Pref="1 bar" Tmax="1075.0" Tmin="250.0">
<floatArray size="7">
50.72389, 6.672267, -2.517167,
10.15934, -0.200675, -427.2115,
130.3973
</floatArray>
</Shomate>
</thermo>
<density units="g/cm3">2.165</density>
</species>
</speciesData>
</ctml>

1
test_problems/cathermo/HMW_graph_VvT/README Normal file
View File

@@ -0,0 +1 @@
Check on the enthalpy routines for HMWSoln

129
test_problems/cathermo/HMW_graph_VvT/TemperatureTable.h Normal file
View File

@@ -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 <vector>
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<double> T;
int numAddedTs;
vector<double> 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

182
test_problems/cathermo/HMW_graph_VvT/V_standalone.cpp Normal file
View File

@@ -0,0 +1,182 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
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);
}

39
test_problems/cathermo/HMW_graph_VvT/output_blessed.txt Normal file
View File

@@ -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

42
test_problems/cathermo/HMW_graph_VvT/runtest Executable file
View File

@@ -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

54
test_problems/cathermo/HMW_graph_VvT/sortAlgorithms.cpp Normal file
View File

@@ -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;
}
}
/*****************************************************/

21
test_problems/cathermo/HMW_graph_VvT/sortAlgorithms.h Normal file
View File

@@ -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