Files
cantera/test/thermo/PengRobinson_Test.cpp
2023-03-06 18:02:07 -06:00

227 lines
7.0 KiB
C++

#include "gtest/gtest.h"
#include "cantera/thermo/PengRobinson.h"
#include "cantera/thermo/ThermoFactory.h"
namespace Cantera
{
class PengRobinson_Test : public testing::Test
{
public:
PengRobinson_Test() {
test_phase = newThermo("../data/thermo-models.yaml", "CO2-PR");
}
//vary the composition of a co2-h2 mixture:
void set_r(const double r) {
vector_fp moleFracs(7);
moleFracs[0] = r;
moleFracs[2] = 1-r;
test_phase->setMoleFractions(&moleFracs[0]);
}
std::shared_ptr<ThermoPhase> test_phase;
};
TEST_F(PengRobinson_Test, chem_potentials)
{
test_phase->setState_TP(298.15, 101325.);
/* Chemical potential should increase with increasing co2 mole fraction:
* mu = mu_0 + RT ln(gamma_k*X_k).
* where gamma_k is the activity coefficient. Run regression test against values
* calculated using the model.
*/
const double expected_result[9] = {
-457338129.70445037,
-457327078.87912911,
-457317214.31077951,
-457308354.65227401,
-457300353.74028891,
-457293092.45485628,
-457286472.73969948,
-457280413.14238912,
-457274845.44186872
};
double xmin = 0.6;
double xmax = 0.9;
int numSteps = 9;
double dx = (xmax-xmin)/(numSteps-1);
vector_fp chemPotentials(7);
for(int i=0; i < numSteps; ++i)
{
set_r(xmin + i*dx);
test_phase->getChemPotentials(&chemPotentials[0]);
EXPECT_NEAR(expected_result[i], chemPotentials[0], 1.e-6);
}
}
TEST_F(PengRobinson_Test, standardConcentrations)
{
EXPECT_DOUBLE_EQ(test_phase->pressure()/(test_phase->temperature()*GasConstant),
test_phase->standardConcentration(0));
EXPECT_DOUBLE_EQ(test_phase->pressure()/(test_phase->temperature()*GasConstant),
test_phase->standardConcentration(1));
}
TEST_F(PengRobinson_Test, setTP)
{
// Check to make sure that the phase diagram is accurately reproduced for a few select isobars
// All sub-cooled liquid:
const double rho1[6] = {
6.6603507723749249e+002,
1.6824762614489907e+002,
1.6248354709581241e+002,
1.5746729362032696e+002,
1.5302217175386241e+002,
1.4902908974486667e+002
};
// Phase change between temperatures 4 & 5:
const double rho2[6] = {
7.5732259810273172e+002,
7.2766981078381912e+002,
6.935475475396446e+002,
6.5227027102964917e+002,
5.9657442842753153e+002,
3.9966973266966875e+002
};
// Supercritical; no discontinuity in rho values:
const double rho3[6] = {
8.0601205067780199e+002,
7.8427655940884574e+002,
7.6105347579146576e+002,
7.3605202492828505e+002,
7.0887891410210011e+002,
6.7898591969734434e+002
};
for(int i=0; i<6; ++i)
{
const double temp = 294 + i*2;
set_r(0.999);
test_phase->setState_TP(temp, 5542027.5);
EXPECT_NEAR(test_phase->density(), rho1[i], 1.e-8);
test_phase->setState_TP(temp, 7388370.);
EXPECT_NEAR(test_phase->density(), rho2[i], 1.e-7);
test_phase->setState_TP(temp, 9236712.5);
EXPECT_NEAR(test_phase->density(), rho3[i], 1.e-8);
}
}
TEST_F(PengRobinson_Test, getPressure)
{
// Check to make sure that the P-R equation is accurately reproduced for a few selected values
/* This test uses CO2 as the only species.
* Values of a_coeff, b_coeff are calculated based on the the critical temperature
* and pressure values of CO2 as follows:
* a_coeff = 0.457235(RT_crit)^2/p_crit
* b_coeff = 0.077796(RT_crit)/p_crit
* The temperature dependent parameter in P-R EoS is calculated as
* \alpha = [1 + \kappa(1 - sqrt{T/T_crit}]^2
* kappa is a function calculated based on the acentric factor.
*/
double a_coeff = 3.958095109E+5;
double b_coeff = 26.62616317/1000;
double acc_factor = 0.228;
double pres_theoretical, kappa, alpha, mv;
const double rho = 1.0737;
//Calculate kappa value
kappa = 0.37464 + 1.54226*acc_factor - 0.26992*acc_factor*acc_factor;
for (int i = 0; i < 15; i++)
{
const double temp = 296 + i * 50;
set_r(1.0);
test_phase->setState_TD(temp, rho);
const double Tcrit = test_phase->critTemperature();
mv = 1 / rho * test_phase->meanMolecularWeight();
//Calculate pressure using Peng-Robinson EoS
alpha = pow(1 + kappa*(1 - sqrt(temp / Tcrit)), 2.0);
pres_theoretical = GasConstant*temp / (mv - b_coeff)
- a_coeff*alpha / (mv*mv + 2*b_coeff*mv - b_coeff*b_coeff);
EXPECT_NEAR(test_phase->pressure(), pres_theoretical, 3);
}
}
TEST_F(PengRobinson_Test, CoolPropValidate)
{
// Validate the P-R EoS in Cantera with P-R EoS from CoolProp
const double rhoCoolProp[10] = {
9.067928191884574,
8.318322900591179,
7.6883521740498155,
7.150504298001246,
6.685330199667018,
6.278630757480957,
5.919763108091383,
5.600572727499541,
5.314694056926007,
5.057077678380463
};
double p = 5e5;
// Calculate density using Peng-Robinson EoS from Cantera
for(int i=0; i<10; i++)
{
const double temp = 300 + i*25;
set_r(1.0);
test_phase->setState_TP(temp, p);
EXPECT_NEAR(test_phase->density(),rhoCoolProp[i],1.e-5);
}
}
TEST(PengRobinson, lookupSpeciesProperties)
{
AnyMap phase_def = AnyMap::fromYamlString(
"{name: test, species: [{gri30.yaml/species: [CO2, CH4, N2]}],"
" thermo: Peng-Robinson}"
);
shared_ptr<ThermoPhase> test(newThermo(phase_def));
// Check for correspondence to values in critical properties "database"
test->setState_TPX(330, 100 * OneAtm, "CH4: 1.0");
EXPECT_NEAR(test->critPressure(), 4.63e+06, 1e-2);
EXPECT_NEAR(test->critTemperature(), 190.7, 1e-3);
test->setState_TPX(330, 180 * OneAtm, "N2: 1.0");
EXPECT_NEAR(test->critPressure(), 3.39e+06, 1e-2);
EXPECT_NEAR(test->critTemperature(), 126.2, 1e-3);
}
TEST(PengRobinson, lookupSpeciesPropertiesMissing)
{
AnyMap phase_def = AnyMap::fromYamlString(
"{name: test, species: [{gri30.yaml/species: [CO2, CH3, N2]}],"
" thermo: Peng-Robinson}"
);
// CH3 is not in the critical properties database, so this should be
// detected as an error
EXPECT_THROW(newThermo(phase_def), CanteraError);
}
TEST(PengRobinson, localCritProperties)
{
// Test calculation based on critical properties stored in the YAML species
// definition, in the "critical-parameters" field
shared_ptr<ThermoPhase> test(newThermo("thermo-models.yaml", "CO2-PR-params"));
test->setState_TPX(400, 1.2e6, "CO2: 1.0");
EXPECT_NEAR(test->critTemperature(), 304.128, 1e-5);
EXPECT_NEAR(test->critPressure(), 7.3773e6, 1e-4);
test->setState_TPX(400, 1.2e6, "H2O: 1.0");
EXPECT_NEAR(test->critTemperature(), 647.096, 1e-5);
EXPECT_NEAR(test->critPressure(), 22.064e6, 1e-4);
}
};