cantera/test_problems/cathermo/HMW_graph_CpvT/Cp_standalone.cpp

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

using namespace std;


/*
 *  Values of A_J/R : tabular form
 *       units  sqrt(kg/gmol), 
 */
double A_JdR(double temp) {
  double retn;
  if (temp == 323.15) {
    retn = 5.50274;
  } else if (temp == 473.15) {
    retn = 24.8263;
  } else {
    printf("A_JdR unknown temp value %g\n", temp);
    exit(-1);
  }
  return retn;
}

double Beta0(double temp, int ifunc) {
  double q1 = 0.0765;
  double q2 = -777.03;
  double q3 = -4.4706;
  double q4 = 0.008946;
  double q5 = -3.3158E-6;
  double retn;
  double tref = 298.15;
  if (ifunc == 0) {
    retn = q1 + q2 * (1.0/temp - 1.0/tref) 
      + q3 * (log(temp/tref)) + q4 * (temp - tref)
      + q5 * (temp * temp - tref * tref);
  } else if (ifunc == 1) {
    retn = (- q2 * 1.0/(temp* temp)  
	    + q3 / temp
	    + q4
	    + 2.0 * temp * q5);
  } else if (ifunc == 2) {
    retn = ( 2.0 * q2 * 1.0/(temp* temp*temp)  
	     - q3 / (temp*temp)
	    + 2.0 * q5);
  } else {
    exit(-1);
  }
  return retn;
}

double Beta1(double temp, int ifunc) {
  double q6 = 0.2664;
  double q9 = 6.1608E-5;
  double q10 = 1.0715E-6;
  double retn;
  double tref = 298.15;
  if (ifunc == 0) {
    retn = q6  + q9 * (temp - tref)
      + q10 * (temp * temp - tref * tref);
  } else if (ifunc == 1) {
    retn = q9 + 2.0 * q10 * temp;
 } else if (ifunc == 2) {
    retn = 2.0 * q10;
  } else {
    exit(-1);
  }
  return retn;
}

double Cphi(double temp, int ifunc) {
  double q11 = 0.00127;
  double q12 = 33.317;
  double q13 = 0.09421;
  double q14 = -4.655E-5;
  double retn;
  double tref = 298.15;
  if (ifunc == 0) {
    retn = q11 + q12 * (1.0/temp - 1.0/tref) 
      + q13 * (log(temp/tref)) + q14 * (temp - tref);
  } else if (ifunc == 1) {
    retn = - q12 / (temp * temp)
      + q13 / temp + q14;
 } else if (ifunc == 2) {
    retn = + 2.0 * q12 / (temp * temp * temp)
      - q13 / (temp * temp) ;
  } else {
    exit(-1);
  }
  return retn;
}

double calc(double temp, double Iionic) {
  /*
   * Gas Constant in J gmol-1 K-1
   */
  double GasConst = 8.314472;

  double Aphi = 0.0;
  if (temp == 323.15) {
    Aphi = 0.4102995331359;
  } else if (temp == 473.15) {
    Aphi = 0.622777;
  } else {
    printf("ERROR: unknown temp\n");
    exit(-1);
  }
  //printf("   Aphi = %g\n", Aphi);

  /*
   * Calculate A_H in J gmol-1 sqrt(kg/gmol)
   */
  double A_J = A_JdR(temp);
  A_J *= GasConst;

  double beta0prime2 = Beta0(temp, 2);
  printf("   beta0prime2 = %g\n", beta0prime2);

  double beta1prime2 = Beta1(temp, 2);
  printf("   beta1prime2 = %g\n", beta1prime2);

  double cphiprime2 = Cphi(temp, 2);
  printf("   Cphiprime2 = %g\n", cphiprime2);

  double beta0prime = Beta0(temp, 1);
  printf("   beta0prime2 = %g\n", beta0prime);

  double beta1prime = Beta1(temp, 1);
  printf("   beta1prime = %g\n", beta1prime);

  double cphiprime = Cphi(temp, 1);
  printf("   Cphiprime = %g\n", cphiprime);

  double vm = 1.0;
  double vx = 1.0;
  double v = vm + vx;
  double m = Iionic;
  double zm = 1.;
  double zx = 1.0;

  double sqrtI = sqrt(Iionic);

  double alpha = 2.0;
  double a2 = alpha * alpha;
  double b = 1.2;

  double Bpmx = beta0prime + 2.0 * beta1prime / (a2* Iionic) *
    (1.0 - (1.0 + alpha * sqrtI) * exp(-alpha*sqrtI) );

  double Bppmx = beta0prime2 + 2.0 * beta1prime2 / (a2* Iionic) *
    (1.0 - (1.0 + alpha * sqrtI) * exp(-alpha*sqrtI) );

  double Cpmx = 0.5 * sqrt(vm * vx) * cphiprime;

  double Cppmx = 0.5 * sqrt(vm * vx) * cphiprime2;

  double Bmx =  Bppmx + 2.0 / temp * Bpmx;
  double Cmx =  Cppmx + 2.0 / temp * Cpmx;

  double phiJ = v * zm * zx * (A_J/(2.*b)) * log(1 + 1.2 * sqrtI) - 
    2 * vm * vx * GasConst * temp * temp * ( m * Bmx + m * m * Cmx);
  phiJ *= 1.0E-3;
  printf("   phiJ = %15.8g kJ/gmolSalt\n", phiJ);

  double molecWeight = 18.01528;

  //double RT = GasConst * temp * 1.0E-3;

  double xo = 1.0 / (molecWeight/1000. * 2 * m + 1.0);
  printf(" no = %g\n", xo);

  return phiJ;
}

int main() {

  printf("Standalone test of the apparent relative molal enthalpy, phiL:\n");
  printf(" (Check against simple formula in Silvester&Pitzer, J. Phys. Chem. 81, 1822 (1977)\n");

  printf("T = 50C\n");
  double Iionic = 6.146;
  printf("Ionic Strength = %g\n", Iionic);

  double res = calc(273.15 + 50., Iionic);
  printf("T = 200C\n");
  printf("Ionic Strength = %g\n", Iionic);

  res = calc(273.15 + 200., Iionic);
  return 0;
}