mirror of
https://github.com/CoolProp/CoolProp.git
synced 2026-01-15 00:48:18 -05:00
111 lines
3.3 KiB
C++
111 lines
3.3 KiB
C++
#include <iostream>
|
|
#include "Eigen/Dense"
|
|
#include "time.h"
|
|
#include "Helmholtz.h"
|
|
#include "CoolProp.h"
|
|
|
|
class EOSFitter;
|
|
#include "Fitter.h"
|
|
#include "DataTypes.h"
|
|
|
|
|
|
int main()
|
|
{
|
|
double n[]={0.0, 0.5586817e-3, 0.4982230e0, 0.2458698e-0, 0.8570145e-3, 0.4788584e-3, -0.1800808e-1, 0.2671641e0, -0.4781652e1, 0.1423987e1, 0.3324062e0, -0.7485907e-2, 0.1017263e-3, -0.5184567e+0, -0.8692288e-1, 0.2057144e+0, -0.5000457e-2, 0.4603262e-3, -0.3497836e-2, 0.6995038e-2, -0.1452184e-1, -0.1285458e-3};
|
|
double d[]={0,2,1,3,6,6,1,1,2,5,2,2,4,1,4,1,2,4,1,5,3,10};
|
|
double t[]={0.0,-1.0/2.0,0.0,0.0,0.0,3.0/2.0,3.0/2.0,2.0,2.0,1.0,3.0,5.0,1.0,5.0,5.0,6.0,10.0,10.0,10.0,18.0,22.0,50.0};
|
|
double c[]={0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,1.0,2.0,2.0,2.0,2.0,2.0,2.0,3.0,3.0,3.0,4.0};
|
|
|
|
std::vector<double> nv(n,n+sizeof(n)/sizeof(double));
|
|
|
|
double mm = Props1SI("R134a","molemass");
|
|
double rhoL, rhoV;
|
|
bool supercritical_T;
|
|
|
|
double Tr = Props1SI("R134a","Treduce");
|
|
EOSFitter* pEOS = new EOSFitterFixedForm(Props1SI("R134a","Treduce"),Props1SI("R134a","rhoreduce")/mm*1000,8.314471);
|
|
EOSFitter &EOS = *pEOS;
|
|
|
|
// ----------------------------
|
|
// Generate "experimental" data
|
|
// ----------------------------
|
|
for (double T = 250; T < 500; T+=10)
|
|
{
|
|
if (T < Tr)
|
|
{
|
|
rhoL = PropsSI("D","T",T,"Q",0,"R134a");
|
|
rhoV = PropsSI("D","T",T,"Q",1,"R134a");
|
|
supercritical_T = false;
|
|
}
|
|
else
|
|
{
|
|
rhoL = -1;
|
|
rhoV = -1;
|
|
supercritical_T = true;
|
|
}
|
|
|
|
for (double rho = 1e-10; rho < 1200; rho *= 1.5)
|
|
{
|
|
if (!supercritical_T && (rho < rhoL && rho > rhoV)){ continue; }
|
|
double p = PropsSI("P","T",T,"D",rho,"R134a");
|
|
double rhobar = rho/mm*1000;
|
|
double cp = PropsSI("C","T",T,"D",rho,"R134a"); // [J/kg/K]; convert to J/mol/K by *mm/1000
|
|
double variance = 1; // TODO; change this
|
|
EOS.linear_data_points.push_back(new PressureDataPoint(pEOS,T,rho/mm*1000,p,variance));
|
|
EOS.nonlinear_data_points.push_back(new SpecificHeatCPDataPoint(pEOS,T,rho/mm*1000,cp*mm/1000,variance*100));
|
|
}
|
|
}
|
|
|
|
// Setup the EOS
|
|
EOS.alphar = phir_power(n,d,t,c,1,21,22);
|
|
|
|
static const double a0[]={
|
|
0.0, //[0]
|
|
-1.019535, //[1]
|
|
9.047135, //[2]
|
|
-1.629789, //[3]
|
|
-9.723916, //[4]
|
|
-3.927170 //[5]
|
|
};
|
|
static const double t0[]={
|
|
0.0, //[0]
|
|
0.0, //[1]
|
|
0.0, //[2]
|
|
0.0, //[3]
|
|
-1.0/2.0, //[4]
|
|
-3.0/4.0 //[5]
|
|
};
|
|
|
|
// phi0=log(delta)+a0[1]+a0[2]*tau+a0[3]*log(tau)+a0[4]*pow(tau,-1.0/2.0)+a0[5]*pow(tau,-3.0/4.0);
|
|
EOS.alpha0.push_back(new phi0_lead(a0[1],a0[2]));
|
|
EOS.alpha0.push_back(new phi0_logtau(a0[3]));
|
|
EOS.alpha0.push_back(new phi0_power(a0,t0,4,5,6));
|
|
|
|
/*for (unsigned int i = 0; i < EOS.nonlinear_data_points.size();i++)
|
|
{
|
|
std::cout << EOS.nonlinear_data_points[i]->residual(nv) << std::endl;
|
|
}*/
|
|
|
|
// Set the coefficients in the preliminary EOS
|
|
EOS.set_n(nv);
|
|
std::cout << format("before fit x2 %g\n",EOS.sum_squares(nv,false));
|
|
// Solve for n without nonlinear terms to get an approximate solution
|
|
EOS.solve_for_n(nv, false);
|
|
std::cout << format("solved for n x2 %g\n",EOS.sum_squares(nv,false));
|
|
EOS.set_n(nv);
|
|
std::cout << format("applied n x2 %g\n",EOS.sum_squares(nv,true));
|
|
|
|
for (int iter = 0; iter < 5; iter++)
|
|
{
|
|
EOS.set_n(nv);
|
|
|
|
// Turn on the nonlinear terms and try again
|
|
EOS.solve_for_n(nv, true);
|
|
|
|
std::cout << nv[1] << " " << nv[2] << std::endl;
|
|
|
|
std::cout << format("iter: %d x2 %g\n",iter, EOS.sum_squares(nv,true));
|
|
}
|
|
|
|
double rr = 0;
|
|
} |