CoolProp
IdealCurves.h
1 #include "AbstractState.h"
2 #include "crossplatform_shared_ptr.h"
3 #include "Solvers.h"
4 #include "CoolPropTools.h"
5 #include <string>
6 
7 namespace CoolProp{
8 
9  class CurveTracer : public FuncWrapper1D
10  {
11  public:
12  AbstractState *AS;
13  double p0, T0, lnT, lnp, rho_guess;
14  std::vector<double> T, p;
15  enum OBJECTIVE_TYPE {OBJECTIVE_INVALID = 0, OBJECTIVE_CIRCLE, OBJECTIVE_T };
16  OBJECTIVE_TYPE obj;
17  CurveTracer(AbstractState *AS, double p0, double T0) : AS(AS), p0(p0), T0(T0), lnT(_HUGE), lnp(_HUGE), rho_guess(_HUGE), obj(OBJECTIVE_INVALID)
18  {
19  this->p.push_back(p0);
20  };
21  void init(){
22  // Solve for Temperature for first point
23  this->obj = OBJECTIVE_T;
24  this->rho_guess = -1;
25  this->T.push_back(Secant(this, T0, 0.001*T0, 1e-10, 100));
26  }
27 
28  virtual double objective(void) = 0;
29 
30  virtual double starting_direction(){
31  return M_PI / 2.0;
32  }
33 
34  double call(double t){
35  if (this->obj == OBJECTIVE_CIRCLE){
36  double T2, P2;
37  this->TPcoords(t, lnT, lnp, T2, P2);
38  this->AS->update(PT_INPUTS, P2, T2);
39  }
40  else{
41  if (this->rho_guess < 0)
42  this->AS->update(PT_INPUTS, this->p[this->p.size()-1], t);
43  else{
44  GuessesStructure guesses;
45  guesses.rhomolar = this->rho_guess;
46  this->AS->update_with_guesses(PT_INPUTS, this->p[this->p.size() - 1], t, guesses);
47  }
48  }
49  double r = this->objective();
50  return r;
51  }
52 
53  void TPcoords(double t, double lnT, double lnp, double &T, double &p){
54  double rlnT = 0.1, rlnp = 0.1;
55  T = exp(lnT + rlnT*cos(t));
56  p = exp(lnp + rlnp*sin(t));
57  }
58 
59  void trace(std::vector<double> &T, std::vector<double> &p)
60  {
61  double t = this->starting_direction();
62  for (int i = 0; i < 1000; ++i){
63  try{
64  this->lnT = log(this->T[this->T.size() - 1]);
65  this->lnp = log(this->p[this->p.size() - 1]);
66  this->obj = OBJECTIVE_CIRCLE;
67  t = Brent(this, t - M_PI / 2.0, t + M_PI / 2.0, DBL_EPSILON, 1e-10, 100);
68  double T2, P2;
69  this->TPcoords(t, this->lnT, this->lnp, T2, P2);
70  this->T.push_back(T2);
71  this->p.push_back(P2);
72  if (this->T[this->T.size() - 1] < this->AS->keyed_output(iT_triple) || this->p[this->p.size() - 1] > 1000 * this->AS->keyed_output(iP_critical)){
73  break;
74  }
75  }
76  catch (std::exception &){
77  break;
78  }
79  }
80  T = this->T;
81  p = this->p;
82  }
83  };
84 
85  class IdealCurveTracer : public CurveTracer{
86  public:
87  IdealCurveTracer(AbstractState *AS, double p0, double T0) : CurveTracer(AS, p0, T0) { init(); };
89  double objective(void){ return this->AS->keyed_output(iZ) - 1; };
90  };
91 
92  class BoyleCurveTracer : public CurveTracer{
93  public:
94  BoyleCurveTracer(AbstractState *AS, double p0, double T0) : CurveTracer(AS, p0, T0) { init(); };
96  double objective(void){
97  double r = (this->AS->p() - this->AS->rhomolar()*this->AS->first_partial_deriv(iP, iDmolar, iT)) / (this->AS->gas_constant()*this->AS->T());
98  return r;
99  };
100  };
102  public:
103  JouleInversionCurveTracer(AbstractState *AS, double p0, double T0) : CurveTracer(AS, p0, T0) { init(); };
105  double objective(void){
106  double r = (this->AS->gas_constant()*this->AS->T() * 1 / this->AS->rhomolar()*this->AS->first_partial_deriv(iP, iT, iDmolar) - this->AS->p()*this->AS->gas_constant() / this->AS->rhomolar()) / POW2(this->AS->gas_constant()*this->AS->T());
107  return r;
108  };
109  };
111  public:
112  JouleThomsonCurveTracer(AbstractState *AS, double p0, double T0) : CurveTracer(AS, p0, T0) { init(); };
114  double objective(void){
115  double dvdT__constp = -this->AS->first_partial_deriv(iDmolar, iT, iP) / POW2(this->AS->rhomolar());
116  double r = this->AS->p() / (this->AS->gas_constant()*POW2(this->AS->T()))*(this->AS->T()*dvdT__constp - 1 / this->AS->rhomolar());
117  return r;
118  };
119  };
120 
121 } /* namespace CoolProp */
double T(void)
Return the temperature in K.
Definition: AbstractState.h:700
Definition: IdealCurves.h:92
double gas_constant(void)
Return the mole-fraction weighted gas constant in J/mol/K.
Definition: AbstractState.cpp:651
double rhomolar(void)
Return the molar density in mol/m^3.
Definition: AbstractState.h:702
CoolPropDbl first_partial_deriv(parameters Of, parameters Wrt, parameters Constant)
The first partial derivative in homogeneous phases.
Definition: AbstractState.h:839
virtual void update_with_guesses(CoolProp::input_pairs input_pair, double Value1, double Value2, const GuessesStructure &guesses)
Update the state using two state variables and providing guess values Some or all of the guesses will...
Definition: AbstractState.h:549
Definition: IdealCurves.h:9
Pressure at the critical point.
Definition: DataStructures.h:62
Definition: IdealCurves.h:85
double objective(void)
Z = 1.
Definition: IdealCurves.h:89
Pressure in Pa, Temperature in K.
Definition: DataStructures.h:227
The mother of all state classes.
Definition: AbstractState.h:70
The compressibility factor Z = p*v/(R*T)
Definition: DataStructures.h:134
double rhomolar
molar density in mol/m^3
Definition: AbstractState.h:34
double objective(void)
dZ/dv|T = 0
Definition: IdealCurves.h:96
double objective(void)
dZ/dT|v = 0
Definition: IdealCurves.h:105
Triple point temperature.
Definition: DataStructures.h:64
double Brent(FuncWrapper1D *f, double a, double b, double macheps, double t, int maxiter)
This function implements a 1-D bounded solver using the algorithm from Brent, R.
Definition: Solvers.cpp:414
Mole-based density.
Definition: DataStructures.h:80
double p(void)
Return the pressure in Pa.
Definition: AbstractState.h:706
Pressure.
Definition: DataStructures.h:74
Definition: IdealCurves.h:110
Definition: Solvers.h:19
This simple class holds the values for guesses for use in some solvers that have the ability to use g...
Definition: AbstractState.h:32
virtual void update(CoolProp::input_pairs input_pair, double Value1, double Value2)=0
Update the state using two state variables.
double objective(void)
dZ/dT|p = 0
Definition: IdealCurves.h:114
This file contains flash routines in which the state is unknown, and a solver of some kind must be us...
Definition: AbstractState.h:19
Temperature.
Definition: DataStructures.h:73
Definition: IdealCurves.h:101
double Secant(FuncWrapper1D *f, double x0, double dx, double ftol, int maxiter)
In the secant function, a 1-D Newton-Raphson solver is implemented.
Definition: Solvers.cpp:286
double keyed_output(parameters key)
Retrieve a value by key.
Definition: AbstractState.cpp:392