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 thermo update

Browse Source
This commit is contained in:
jeremy karpowski 2025-02-03 09:06:30 +01:00
parent e36c9b9bcf
commit 0710fef5f8
5 changed files with 115 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
include/cantera/thermo/IdealGasPhase.h
View File

@ -366,7 +366,7 @@ public:
* @param p Pressure (Pa)
* @since New in %Cantera 3.0.
*/
void setState_DP(double rho, double p) override {
void setState_DP(double rho, double p, double Tguess=300) override {
if (p <= 0) {
throw CanteraError("IdealGasPhase::setState_DP",
"pressure must be positive");

55
include/cantera/thermo/PengRobinson.h
View File

@ -203,6 +203,7 @@ protected:
* These are stored internally.
*/
double daAlpha_dT() const;
double daAlpha_dT(double& temp) const;
//! Calculate second derivative @f$ d^2(a \alpha)/dT^2 @f$
/*!
@ -211,6 +212,58 @@ protected:
double d2aAlpha_dT2() const;
public:
void setState_DP(double rho, double p, double Tguess=300) override {
if (p <= 0) {
throw CanteraError("IdealGasPhase::setState_DP",
"pressure must be positive");
}
//"Yeah its this one"
setDensity(rho);
setForcedSolutionBranch(-2);
double mV = meanMolecularWeight() / rho;
// Update individual alpha
double aCalc, bCalc, aAlpha;
int iteration = 0;
double T_new = Tguess; //(p + term2) * denum / GasConstant;
double p_iter=p;
while (iteration < 300) {
// Update a, b, and alpha based on the current guess for T
updateMixingExpressions(T_new, aCalc, bCalc, aAlpha);
// Calculate p using PR EoS
double denom = mV * mV + 2 * mV * bCalc - bCalc * bCalc;
double mV_b = mV - bCalc;
double term2 = (aAlpha) / denom;
p_iter = (GasConstant*T_new)/mV_b - term2;
//std::cout << "T_ : "<<T_new << "p : "<<p <<std::endl;
// Check if the result is within the specified tolerance
double error =(p_iter-p)/p;
//std::cout << "error: "<< error <<" dpdT: "<<m_dpdT/p <<" T: "<<T_new<<" p: "<<p << " p_iter: "<<p_iter<<" rho: "<<rho<<std::endl;
//
if (fabs(error)<1e-10) {
break;
}
m_dpdT = (GasConstant / mV_b - daAlpha_dT(T_new) / denom);
T_new = T_new - (error/(m_dpdT/p));
iteration++;
}
if (iteration==300)
{
throw CanteraError("PengRobinson::DP did not converge at",
"negative temperature T = {} p= {} rho={} X1={} X2={}", T_new, p, rho,moleFractions_[0],moleFractions_[1]);
}
// if ((rho > 50) && (rho<1000))
// {
// std::cout <<"T_new: "<<T_new <<" iteration "<< iteration<<" forcedState "<< forcedState_<<std::endl;
// }
setDensity(rho);
setTemperature(T_new);
setState_TP(T_new,p);
updateMixingExpressions();
}
double isothermalCompressibility() const override;
double thermalExpansionCoeff() const override;
@ -229,7 +282,7 @@ public:
* the internal numbers based on the state of the object.
*/
void updateMixingExpressions() override;
void updateMixingExpressions(double& temp,double& aCalc, double& bCalc, double& aAlphaCalc);
//! Calculate the @f$ a @f$, @f$ b @f$, and @f$ \alpha @f$ parameters given the temperature
/*!
* This function doesn't change the internal state of the object, so it is a

2
include/cantera/thermo/ThermoPhase.h
View File

@ -1334,7 +1334,7 @@ public:
* @param p Pressure (Pa)
* @since New in %Cantera 3.0.
*/
virtual void setState_DP(double rho, double p) {
virtual void setState_DP(double rho, double p, double Tguess=300) {
throw NotImplementedError("ThermoPhase::setState_DP");
}

17
src/thermo/MixtureFugacityTP.cpp
View File

@ -835,8 +835,14 @@ int MixtureFugacityTP::solveCubic(double T, double pres, double a, double b,
//check if y = h
if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) {
if (disc > 1e-10) {
std::cout<<"yn: "<<yn<<" h: "<<h<<" disc: "<<disc<<std::endl;
writelog("MixtureFugacityTP::solveCubic(T ={}, p ={}):"
"X1: {}, X2: {}\n",
T, pres
, moleFractions_[0], moleFractions_[1]);
throw CanteraError("MixtureFugacityTP::solveCubic",
"value of yN and h are too high, unrealistic roots may be obtained");
"value of yN and h are too high, unrealistic roots may be obtained"
);
}
disc = 0.0;
}
@ -893,8 +899,9 @@ int MixtureFugacityTP::solveCubic(double T, double pres, double a, double b,
for (int j = 0; j < 3; j++) {
if (j != i && fabs(Vroot[i] - Vroot[j]) < 1.0E-4 * (fabs(Vroot[i]) + fabs(Vroot[j]))) {
writelog("MixtureFugacityTP::solveCubic(T ={}, p ={}):"
" WARNING roots have merged: {}, {}\n",
T, pres, Vroot[i], Vroot[j]);
" WARNING roots have merged: {}, {}\n"
"X1: {}, X2: {}\n",
T, pres, Vroot[i], Vroot[j], moleFractions_[0], moleFractions_[1]);
}
}
}
@ -955,7 +962,9 @@ int MixtureFugacityTP::solveCubic(double T, double pres, double a, double b,
}
if ((fabs(res) > 1.0E-14) && (fabs(res) > 1.0E-14 * fabs(dresdV) * fabs(Vroot[i]))) {
writelog("MixtureFugacityTP::solveCubic(T = {}, p = {}): "
"WARNING root didn't converge V = {}", T, pres, Vroot[i]);
"WARNING root didn't converge V = {}\n"
"X1: {}, X2: {}\n",
T, pres, Vroot[i],moleFractions_[0],moleFractions_[1]);
writelogendl();
}
}

46
src/thermo/PengRobinson.cpp
View File

@ -673,6 +673,28 @@ void PengRobinson::updateMixingExpressions()
calculateAB(m_a,m_b,m_aAlpha_mix);
}
void PengRobinson::updateMixingExpressions(double& temp,double& aCalc, double& bCalc, double& aAlphaCalc)
{
//double temp = temperature();
// Update individual alpha
for (size_t j = 0; j < m_kk; j++) {
double critTemp_j = speciesCritTemperature(m_a_coeffs(j,j), m_b_coeffs[j]);
double sqt_alpha = 1 + m_kappa[j] * (1 - sqrt(temp / critTemp_j));
m_alpha[j] = sqt_alpha*sqt_alpha;
}
// Update aAlpha_i, j
for (size_t i = 0; i < m_kk; i++) {
for (size_t j = 0; j < m_kk; j++) {
m_aAlpha_binary(i, j) = sqrt(m_alpha[i] * m_alpha[j]) * m_a_coeffs(i,j);
}
}
calculateAB(aCalc,bCalc,aAlphaCalc);
}
void PengRobinson::calculateAB(double& aCalc, double& bCalc, double& aAlphaCalc) const
{
bCalc = 0.0;
@ -710,6 +732,30 @@ double PengRobinson::daAlpha_dT() const
return daAlphadT;
}
double PengRobinson::daAlpha_dT(double& temp) const
{
double daAlphadT = 0.0, k, Tc, sqtTr, coeff1, coeff2;
for (size_t i = 0; i < m_kk; i++) {
// Calculate first derivative of alpha for individual species
Tc = speciesCritTemperature(m_a_coeffs(i,i), m_b_coeffs[i]);
sqtTr = sqrt(temp / Tc);
coeff1 = 1 / (Tc*sqtTr);
coeff2 = sqtTr - 1;
k = m_kappa[i];
m_dalphadT[i] = coeff1 * (k*k*coeff2 - k);
}
// Calculate mixture derivative
for (size_t i = 0; i < m_kk; i++) {
for (size_t j = 0; j < m_kk; j++) {
daAlphadT += moleFractions_[i] * moleFractions_[j] * 0.5
* m_aAlpha_binary(i, j)
* (m_dalphadT[i] / m_alpha[i] + m_dalphadT[j] / m_alpha[j]);
}
}
return daAlphadT;
}
double PengRobinson::d2aAlpha_dT2() const
{
for (size_t i = 0; i < m_kk; i++) {