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 ForceState option to thermo-phase

Browse Source
This commit is contained in:
jeremy karpowski 2025-02-03 10:56:23 +01:00
parent 94cf35e3a4
commit fa6d4d9734
4 changed files with 63 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -88,7 +88,7 @@ public:
* -1 means gas. The value -2 means unrestricted. Values of zero or
* greater refer to species dominated condensed phases.
*/
void setForcedSolutionBranch(int solnBranch);
void setForcedSolutionBranch(int solnBranch) override;
//! Report the solution branch which the solution is restricted to
/*!

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

@ -212,58 +212,7 @@ 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();
}
void setState_DP(double rho, double p, double Tguess=300) override;
double isothermalCompressibility() const override;
double thermalExpansionCoeff() const override;

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

@ -1405,6 +1405,14 @@ public:
throw NotImplementedError("ThermoPhase::setState_DP",
"Not implemented for phase type '{}'", type());
}
//! Set the solution branch to force the ThermoPhase to exist on one branch
//! or another
/*!
* @param solnBranch Branch that the solution is restricted to. the value
* -1 means gas. The value -2 means unrestricted. Values of zero or
* greater refer to species dominated condensed phases.
*/
virtual void setForcedSolutionBranch(int solnBranch){};
//! Set the state using an AnyMap containing any combination of properties
//! supported by the thermodynamic model

53
src/thermo/PengRobinson.cpp
View File

@ -838,4 +838,57 @@ AnyMap PengRobinson::getAuxiliaryData()
return parameters;
}
void PengRobinson::setState_DP(double rho, double p, double Tguess)
{
if (p <= 0 || rho<=0 ) {
throw CanteraError("PengRobinson::setState_DP",
"pressure and density must be positive");
}
//"Yeah its this one"
setDensity(rho);
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();
}
}