CoolProp  4.2.5
An open-source fluid property and humid air property database
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AceticAcid.cpp
Go to the documentation of this file.
1 #include "CoolProp.h"
2 #include <vector>
3 #include "CPExceptions.h"
4 #include "FluidClass.h"
5 #include "AceticAcid.h"
6 #include "Solvers.h"
7 
9 class DensityTpResids : public FuncWrapper1D
10 {
11 private:
12  double p,T;
13  Fluid *pFluid;
14 public:
15  DensityTpResids(Fluid *pFluid, double T, double p){this->pFluid = pFluid; this->p = p; this->T = T;};
16  ~DensityTpResids(){};
17 
18  double call(double rho)
19  {
20  return this->p - pFluid->pressure_Trho(T,rho);
21  }
22 };
23 
24 AceticAcidClass::AceticAcidClass()
25 {
26  double n[] = {0, -0.15624834164583e1, -0.874703669570960e0, 0.46968858010355e1, 0.97367136204905e-2, -0.49055972708048e-2, 0.24499997808125e2, -0.31443235067567e2, -0.13768156877983e1, 0.14849435860881e1, 0.11374909453775e1, -0.26039791873344e1, -0.30484923493199e-1, 0.53316386834696e1, -0.56733952193640e1, -0.126785566440530e0};
27  double d[] = {0, 1, 1, 2, 2, 6, 3, 3, 3, 4, 4, 5, 5, 5, 5, 2};
28  double t[] = {0, -1.000, 1.375, 1.000, 1.375, 0.750, -0.250, 0.000, 2.250, 0.125, 2.125, 1.250, 2.250, 2.125, 2.375, 14.000};
29  double l[] = {0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3};
30 
31  // Critical parameters
32  crit.rho = 351;
33  crit.T = 590.70;
34  crit.v = 1.0/crit.rho;
35 
36  // Other fluid parameters
37  params.molemass = 60.05196;
38  params.Ttriple = 289.8;
39  params.ptriple = _HUGE;
40  params.accentricfactor = 0.48092;
41  params.R_u = 8.314472;
42 
43  // Limits of EOS
44  limits.Tmin = params.Ttriple;
45  limits.Tmax = 500.0;
46  limits.pmax = 100000.0;
47  limits.rhomax = 1000000.0*params.molemass;
48 
49  // Residual part
50  phirlist.push_back(new phir_power(n,d,t,l,1,15,16));
51  double m = 1.01871348, vbarn = 0.444215309e-1,kappabar = 0.109117041e-4 , epsilonbar = 12.2735737;
52  phirlist.push_back(new phir_SAFT_associating_1(m, epsilonbar, vbarn, kappabar));
53 
54  // Ideal-gas part
55  phi0list.push_back(new phi0_lead(-3.94616949, 5.48487930));
56  phi0list.push_back(new phi0_logtau(3.66766530));
57  phi0list.push_back(new phi0_power(-0.210687796, -1));
58  phi0list.push_back(new phi0_power(-0.781330239, -2));
59  phi0list.push_back(new phi0_power(0.130979005, -3));
60  phi0list.push_back(new phi0_Planck_Einstein( 6.28891793, 2.09502491));
61 
62  // Set the critical pressure based on the evaluation of the EOS
63  reduce = crit;
64  crit.p = PressureUnit(pressure_Trho(crit.T,crit.rho+1e-10), UNIT_PA);
65 
66  name.assign("AceticAcid");
67  aliases.push_back(std::string("ACETICACID"));
68  REFPROPname.assign("N/A");
69 
70  BibTeXKeys.EOS = "Piazza-FPE-2011";
71  BibTeXKeys.SURFACE_TENSION = "Mulero-JPCRD-2012";
72 
73  params.CAS = "64-19-7";
74 
75  //double w = 6.67228479e-09*Tc*Tc*Tc-7.20464352e-06*Tc*Tc+3.16947758e-03*Tc-2.88760012e-01;
76  //double q = -6.08930221451*w-5.42477887222;
77  //double pt = exp(q*(Tc/Tt-1))*pc;
78 
79  /*if (1)
80  {
81  double rhoL,rhoV;
82 
83  std::string errstr;
84  double T = params.Ttriple;
85  double p0 = pt;
86  double gibbsL,gibbsV;
87  double rhoLstore = _HUGE;
88  for (double rho = crit.rho; rho < crit.rho*7; rho += crit.rho/5)
89  {
90  DensityTpResids DTPR = DensityTpResids(this,T,p0);
91  std::string errstr;
92  try{
93  rhoL = density_Tp(T,p0,rho);
94  gibbsL = gibbs_Trho(T,rhoL);
95  double pp = pressure_Trho(T,rhoL);
96  if (dpdrho_Trho(T,rhoL)>0){
97  rhoLstore = rhoL;
98  };
99  double rr = 0;
100  }
101  catch (std::exception &)
102  {
103  }
104  }
105  double rhoVstore = _HUGE;
106  for (double rho = crit.rho; rho > 1e-16; rho /= 1.5)
107  {
108  try{
109  rhoV = density_Tp(T,p0,rho);
110  double pp = pressure_Trho(T,rhoV);
111  gibbsV = gibbs_Trho(T,rhoV);
112  if (dpdrho_Trho(T,rhoV)>0){
113  rhoVstore = rhoV;
114  };
115  double rr = 0;
116  }
117  catch (std::exception &)
118  {
119  }
120  }
121  double p;
122  rhoL = rhoLstore; rhoV = rhoVstore;
123  FILE *fp;
124  fp = fopen("aceticancillary.csv","w");
125  fclose(fp);
126  while (T < Tc)
127  {
128  if (T>486)
129  {
130  rhosatPure_Akasaka(T, &rhoL, &rhoV, &p, 0.1, true);
131  }
132  else
133  {
134  rhosatPure(T, &rhoL, &rhoV, &p, 1.0, true);
135  }
136  fp = fopen("aceticancillary.csv","a+");
137  fprintf(fp,"%g,%g,%g,%g\n", T, p, rhoL, rhoV);
138  fclose(fp);
139  if (T < 580)
140  {
141  T += 1;
142  }
143  else
144  {
145  T += 0.001;
146  }
147  }
148  fclose(fp);
149  double rr = 0;
150  }*/
151 }
152 
153 double AceticAcidClass::psat(double T)
154 {
155  // Max error is 0.104268020681 % between 290.0 and 588.0 K
156  const double t[] = {0, 0.065, 0.367, 1.0, 4.0, 4.166666666666667};
157  const double N[] = {0, 0.03788266794177342, -0.3014667718204638, -7.465403420229625, 51.68803799661326, -61.422675342107055};
158  double summer = 0, theta;
159  theta = 1 - T/crit.T;
160  for (int i=1; i <= 5; i++)
161  {
162  summer += N[i]*pow(theta, t[i]);
163  }
164  return 5780000*exp(crit.T/T*summer);
165 }
166 
167 double AceticAcidClass::rhosatL(double T)
168 {
169  // Max error is 0.24407295586 % between 290.0 and 590.0 K
170  const double t[] = {0, 0.13, 0.3585, 0.361, 0.365, 0.8333333333333334};
171  const double N[] = {0, -8.839463207863323, 73076.09152165553, -118684.81658933076, 45623.75583577903, -3.5208689421322545};
172  double summer=0,theta;
173  theta = 1 - T/crit.T;
174  for (int i=1; i <= 5; i++)
175  {
176  summer += N[i]*pow(theta,t[i]);
177  }
178  double val = crit.rho*(1+summer);
179  return val;
180 }
181 
182 double AceticAcidClass::rhosatV(double T)
183 {
184  // Max error is 1.16918454075 % between 290.0 and 590.0 K
185  const double t[] = {0, 0.14600000000000002, 0.38249999999999995, 0.39349999999999996, 0.39749999999999996, 5.333333333333333};
186  const double N[] = {0, 7.738193007082631, -6937.911692305578, 25918.893876023878, -18994.78347845314, -14.42897045481017};
187  double summer=0,theta;
188  theta=1-T/reduce.T;
189  for (int i=1; i<=5; i++)
190  {
191  summer += N[i]*pow(theta,t[i]);
192  }
193  double val = crit.rho*exp(crit.T/T*summer);
194  return val;
195 }
196 
197 #ifndef DISABLE_CATCH
198 #include "Catch/catch.hpp"
199 TEST_CASE("Acetic acid validation from Piazza, FPE, 2011","[aceticacid],[validation]")
200 {
201  // Temperature (K) Density (kg/m3) Pressure (MPa) Compress. factor (Z) Enthalpy (kJ/kg) Entropy (kJ/kgK) cv (kJ/kgK) cp (kJ/kgK) Sound speed (m/s)
202  double data[10][9] = {
203  {290,1e-16,1e-16,1.00000,669.46135,0,0.89579,1.03425,215.30848},
204  {290,0.010,0.23586e-3,0.58743,238.11539,1.05462,4.64706,5.35799,158.69275},
205  {290,0.025,0.55619e-3,0.55408,203.44932,0.86742,3.46806,3.91294,154.08504},
206  {290,1060.,9.5458,0.22429,-218.59529,-0.66978,1.66179,1.99658,1221.38239},
207  {290,1070.,22.580,0.52557,-209.68922,-0.68127,1.65651,1.98783,1278.95824},
208  {450,1e-16,1e-16,1.00000,868.66312,0,1.31313,1.45159,262.43829},
209  {450,1.0,0.55569e-1,0.89189,760.34138,1.94165,3.74492,4.41190,244.91797},
210  {450,6.0,0.26514,0.70925,599.36342,1.40868,4.43708,5.53182,212.41678},
211  {450,870.,2.7252,0.05028,152.00058,0.35592,2.48256,2.74525,536.47030},
212  {450,880.,5.6865,0.10372,153.23341,0.35114,2.46536,2.73491,607.71381}};
213 
214  int N = 10;
215 
216  SECTION("P")
217  {
218  for (int i = 0; i < N; i++)
219  {
220  double T = data[i][0], rho = data[i][1], p = data[i][2]*1e6;
221  double PCP = PropsSI("P","T",T,"D",rho,"AceticAcid");
222  if (rho < 1e-14){p = 1; PCP = 1;}
223  CAPTURE(rho);
224  CAPTURE(p);
225  CAPTURE(PCP);
226  CHECK(fabs(PCP/p-1) < 1e-4);
227  }
228  }
229  SECTION("H")
230  {
231  for (int i = 0; i < N; i++)
232  {
233  double T = data[i][0], rho = data[i][1], h = data[i][4]*1e3;
234  double HCP = PropsSI("H","T",T,"D",rho,"AceticAcid");
235  CAPTURE(h);
236  CAPTURE(HCP);
237  CHECK(fabs(HCP-h) < 1);
238  }
239  }
240  SECTION("S")
241  {
242  for (int i = 0; i < N; i++)
243  {
244  double T = data[i][0], rho = data[i][1], s = data[i][5]*1e3;
245  double SCP = PropsSI("S","T",T,"D",rho,"AceticAcid");
246  if (rho < 1e-14){s = 1; SCP = 1;}
247  CAPTURE(s);
248  CAPTURE(SCP);
249  CHECK(fabs(SCP-s) < 0.01);
250  }
251  }
252  SECTION("O")
253  {
254  for (int i = 0; i < N; i++)
255  {
256  double T = data[i][0], rho = data[i][1], cv = data[i][6]*1e3;
257  double OCP = PropsSI("O","T",T,"D",rho,"AceticAcid");
258  CAPTURE(cv);
259  CAPTURE(OCP);
260  CHECK(fabs(OCP-cv) < 1e-2);
261  }
262  }
263  SECTION("C")
264  {
265  for (int i = 0; i < N; i++)
266  {
267  double T = data[i][0], rho = data[i][1], cp = data[i][7]*1e3;
268  double CCP = PropsSI("C","T",T,"D",rho,"AceticAcid");
269  CAPTURE(cp);
270  CAPTURE(CCP);
271  CHECK(fabs(CCP-cp) < 1e-2);
272  }
273  }
274  SECTION("A")
275  {
276  for (int i = 0; i < N; i++)
277  {
278  double T = data[i][0], rho = data[i][1], a = data[i][8];
279  double ACP = PropsSI("A","T",T,"D",rho,"AceticAcid");
280  CAPTURE(a);
281  CAPTURE(ACP);
282  CHECK(fabs(ACP-a) < 1e-3);
283  }
284  }
285 }
286 TEST_CASE("Check acetic derivatives residual","[aceticacid],[helmholtz]")
287 {
288  AceticAcidClass Acetic = AceticAcidClass();
289  double eps = sqrt(DBL_EPSILON);
290  SECTION("dDelta")
291  {
292  double ANA = Acetic.dphir_dDelta(0.5, 0.5);
293  double NUM = (Acetic.phir(0.5, 0.5+eps) - Acetic.phir(0.5,0.5-eps))/(2*eps);
294  REQUIRE(fabs(NUM-ANA) < 1e-6);
295  }
296  SECTION("dTau")
297  {
298  double ANA = Acetic.dphir_dTau(0.5, 0.5);
299  double NUM = (Acetic.phir(0.5+eps, 0.5) - Acetic.phir(0.5-eps,0.5))/(2*eps);
300  REQUIRE(fabs(NUM-ANA) < 1e-6);
301  }
302  SECTION("dDelta2")
303  {
304  double ANA = Acetic.d2phir_dDelta2(0.5, 0.5);
305  double NUM = (Acetic.dphir_dDelta(0.5, 0.5+eps) - Acetic.dphir_dDelta(0.5,0.5-eps))/(2*eps);
306  REQUIRE(fabs(NUM-ANA) < 1e-6);
307  }
308  SECTION("dTau2")
309  {
310  double ANA = Acetic.d2phir_dTau2(0.5, 0.5);
311  double NUM = (Acetic.dphir_dTau(0.5+eps, 0.5) -Acetic.dphir_dTau(0.5-eps,0.5))/(2*eps);
312  REQUIRE(fabs(NUM-ANA) < 1e-6);
313  }
314  SECTION("dDeltadTau")
315  {
316  double ANA = Acetic.d2phir_dDelta_dTau(0.5, 0.5);
317  double NUM = (Acetic.dphir_dTau(0.5, 0.5+eps) - Acetic.dphir_dTau(0.5,0.5-eps))/(2*eps);
318  REQUIRE(fabs(NUM-ANA) < 1e-6);
319  }
320 }
321 
322 TEST_CASE("Check acetic derivatives ideal","[aceticacid],[helmholtz]")
323 {
324  AceticAcidClass Acetic = AceticAcidClass();
325  double eps = sqrt(DBL_EPSILON);
326 
327  SECTION("dDelta")
328  {
329  double ANA = Acetic.dphi0_dDelta(0.5, 0.5);
330  double NUM = (Acetic.phi0(0.5, 0.5+eps) - Acetic.phi0(0.5,0.5-eps))/(2*eps);
331  REQUIRE(fabs(NUM-ANA) < 1e-6);
332  }
333  SECTION("dTau")
334  {
335  double ANA = Acetic.dphi0_dTau(0.5, 0.5);
336  double NUM = (Acetic.phi0(0.5+eps, 0.5) - Acetic.phi0(0.5-eps,0.5))/(2*eps);
337  REQUIRE(fabs(NUM-ANA) < 1e-6);
338  }
339  SECTION("dDelta2")
340  {
341  double ANA = Acetic.d2phi0_dDelta2(0.5, 0.5);
342  double NUM = (Acetic.dphi0_dDelta(0.5, 0.5+eps) - Acetic.dphi0_dDelta(0.5,0.5-eps))/(2*eps);
343  REQUIRE(fabs(NUM-ANA) < 1e-6);
344  }
345  SECTION("dTau2")
346  {
347  double ANA = Acetic.d2phi0_dTau2(0.5, 0.5);
348  double NUM = (Acetic.dphi0_dTau(0.5+eps, 0.5) -Acetic.dphi0_dTau(0.5-eps,0.5))/(2*eps);
349  REQUIRE(fabs(NUM-ANA) < 1e-6);
350  }
351  SECTION("dDeltadTau")
352  {
353  double ANA = Acetic.d2phi0_dDelta_dTau(0.5, 0.5);
354  double NUM = (Acetic.dphi0_dTau(0.5, 0.5+eps) - Acetic.dphi0_dTau(0.5,0.5-eps))/(2*eps);
355  REQUIRE(fabs(NUM-ANA) < 1e-6);
356  }
357 }
358 
359 #endif
DensityTpResids::DensityTpResids
DensityTpResids(Fluid *pFluid, double T, double p)
Definition: AceticAcid.cpp:15
AceticAcidClass::rhosatL
double rhosatL(double)
Definition: AceticAcid.cpp:167
CoolProp.h
Fluid::phirlist
std::vector< phi_BC * > phirlist
Definition: FluidClass.h:178
phi0_Planck_Einstein
Definition: Helmholtz.h:600
Fluid::d2phi0_dTau2
virtual double d2phi0_dTau2(double tau, double delta)
Definition: FluidClass.cpp:411
Fluid::d2phir_dDelta_dTau
virtual double d2phir_dDelta_dTau(double tau, double delta)
Definition: FluidClass.cpp:329
Fluid::limits
struct FluidLimits limits
Definition: FluidClass.h:219
BibTeXKeysStruct::EOS
std::string EOS
Definition: FluidClass.h:120
phir_power
Definition: Helmholtz.h:91
AceticAcid.h
CriticalStruct::p
PressureUnit p
Definition: FluidClass.h:50
UNIT_PA
Definition: GlobalConstants.h:36
Fluid::name
std::string name
A container to hold the cache for residual Helmholtz derivatives.
Definition: FluidClass.h:151
Fluid::phi0
virtual double phi0(double tau, double delta)
Definition: FluidClass.cpp:381
Solvers.h
TEST_CASE
TEST_CASE("Acetic acid validation from Piazza, FPE, 2011","[aceticacid],[validation]")
Definition: AceticAcid.cpp:199
Fluid::pressure_Trho
double pressure_Trho(double T, double rho)
Definition: FluidClass.cpp:449
DensityTpResids
A stub class to do the density(T,p) calculations for near the critical point using Brent solver...
Definition: FluidClass.cpp:129
Fluid::aliases
std::vector< std::string > aliases
The REFPROP-compliant name if REFPROP-"name" is not a compatible fluid name. If not included...
Definition: FluidClass.h:153
Fluid::dphi0_dTau
virtual double dphi0_dTau(double tau, double delta)
Definition: FluidClass.cpp:403
phi0_logtau
Definition: Helmholtz.h:568
DBL_EPSILON
#define DBL_EPSILON
Definition: Helmholtz.cpp:10
Fluid::reduce
struct CriticalStruct reduce
A pointer to the point that is used to reduce the T and rho for EOS.
Definition: FluidClass.h:222
FluidLimits::Tmax
double Tmax
Definition: FluidClass.h:54
Fluid
Fluid is the abstract base class that is employed by all the other fluids.
Definition: FluidClass.h:147
CAPTURE
#define CAPTURE(msg)
Definition: catch.hpp:8074
Fluid::dphir_dTau
virtual double dphir_dTau(double tau, double delta)
Definition: FluidClass.cpp:295
PropsSI
double PropsSI(std::string Output, std::string Name1, double Prop1, std::string Name2, double Prop2, std::string Ref)
Definition: CoolProp.cpp:886
Fluid::dphi0_dDelta
virtual double dphi0_dDelta(double tau, double delta)
Definition: FluidClass.cpp:389
Fluid::REFPROPname
std::string REFPROPname
The name of the fluid.
Definition: FluidClass.h:152
DensityTpResids::call
double call(double rho)
Definition: AceticAcid.cpp:18
BibTeXKeysStruct::SURFACE_TENSION
std::string SURFACE_TENSION
Definition: FluidClass.h:126
phi0_lead
Definition: Helmholtz.h:493
phi0_power
Definition: Helmholtz.h:696
FluidLimits::pmax
double pmax
Definition: FluidClass.h:54
params
params
Definition: GlobalConstants.h:11
Fluid::d2phir_dDelta2
virtual double d2phir_dDelta2(double tau, double delta)
Definition: FluidClass.cpp:276
Fluid::crit
struct CriticalStruct crit
Definition: FluidClass.h:218
AceticAcidClass::rhosatV
double rhosatV(double)
Definition: AceticAcid.cpp:182
Fluid::BibTeXKeys
BibTeXKeysStruct BibTeXKeys
Definition: FluidClass.h:175
Fluid::phir
virtual double phir(double tau, double delta)
Definition: FluidClass.cpp:237
REQUIRE
REQUIRE(fabs(CPWater.p()-RPWater.p())< 1e-4)
Fluid::d2phi0_dDelta2
virtual double d2phi0_dDelta2(double tau, double delta)
Definition: FluidClass.cpp:396
Fluid::phi0list
std::vector< phi_BC * > phi0list
A vector of instances of the phi_BC classes for the residual Helmholtz energy contribution.
Definition: FluidClass.h:179
CPExceptions.h
Fluid::d2phi0_dDelta_dTau
virtual double d2phi0_dDelta_dTau(double tau, double delta)
Definition: FluidClass.cpp:427
Fluid::d2phir_dTau2
virtual double d2phir_dTau2(double tau, double delta)
Definition: FluidClass.cpp:312
SECTION
#define SECTION(name, description)
Definition: catch.hpp:8088
Fluid::dphir_dDelta
virtual double dphir_dDelta(double tau, double delta)
Definition: FluidClass.cpp:257
AceticAcidClass::psat
double psat(double)
Definition: AceticAcid.cpp:153
FluidLimits::Tmin
double Tmin
Definition: FluidClass.h:54
FluidClass.h
FluidLimits::rhomax
double rhomax
Definition: FluidClass.h:54
CHECK
#define CHECK(expr)
Definition: catch.hpp:8058
phir_SAFT_associating_1
Definition: Helmholtz.h:474
Generated on Sun Aug 10 2014 16:02:40 for CoolProp by   doxygen 1.8.7