Calculating the liquid density at saturation

Posted on February 6, 2019 by paolog

This is the full source code required to calculate with LIBPF the liquid density at saturation for an equimolar mixture of ethane and propane at pressures between 1 and 30 atm:

/** @file Qsatliqrho.cc
 @brief Calculate the density of a saturated liquid
 
 This file is part of LIBPF
 All rights reserved; do not distribute without permission.
 @author (C) Copyright 2019 Paolo Greppi simevo s.r.l.
 */

/* SYSTEM/GENERAL PURPOSE INCLUDES */
//
#include 

/* LIBPF INCLUDES */
//
#include 
#include 
#include  // for initializeLibpf and uninitializeLibpf
#include 
#include 
#include 
#include  // for components and NCOMPONENTS
#include 
#include 
#include  // register all stream types with the model factory
#include 

/* LOCAL INCLUDES */
//

/* FORWARD REFERENCES */
//

/* MACROS */
// avoid macros!!
//

/* GLOBAL CONSTANTS */
//
static const int verbosityFile = 0;

/* GLOBAL VARIABLES */
//

/*====================================================================================================================*/
/*==  main  ==========================================================================================================*/
/*====================================================================================================================*/

int main(int /* argc */, char * /*argv */ []) {
  const int verbosityLocal = 0;
  const int verbosityInstance = 0;
  diagnostic(0, "Entered");
  try {
    // initializations
    initializeLibpf();

    components.addcomp(new purecomps::ethane);
    components.addcomp(new purecomps::propane);

    // Pcsaft i, m, sigma, epsilon
    Phase::PcsaftParameters::resize(2);
    // ethane
    Phase::PcsaftParameters::setparameters(0, 1.6069, 3.5206, 191.42);
    // propane
    Phase::PcsaftParameters::setparameters(1, 2.0020, 3.6184, 208.11);
    
    NodeFactory nodeFactory;
    // try with: StreamPcsaftLiquidVapor, StreamKwongRedlichSoaveLiquidVapor, StreamPengRobinsonLiquidVapor
    Stream *feed = my_cast(nodeFactory.create("StreamPengRobinsonLiquidVapor", Libpf::Persistency::Defaults("test")), CURRENT_FUNCTION);
    feed->Tphase->Q("x[C2H6]").set(0.5);
    feed->Tphase->Q("x[C3H8]").set(0.5);
    feed->flashoption = "PA";
    feed->flowoption= "Nx";
    feed->Q("Vphase.fraction").set(Zero);

    for (Value P(1.0, "atm"); P < Value(30.0, "atm"); P += Value(0.1, "atm")) {
      feed->P.set(P);
      feed->resetErrors();
      feed->calculate();
      Value rhol(feed->Q("Lphase.rho"));
      diagnostic(0, "P = " << P << " Tb = " << feed->T << " rhol = " << rhol << " errors = " << feed->errors.size())
    }
    delete feed;

    uninitializeLibpf();
    return 0;
  } // try
  catch (Error &e) {
    uninitializeLibpf();
    e.append(CURRENT_FUNCTION);
    std::cerr << "****************************** Fatal LIBPF error! ******************************" << std::endl;
    std::cerr << e.message() << std::endl;
    return -1;
  } // catch
} // main

Results with Redlich-Kwong-Soave equation of state: