bwrs.hpp

#pragma once
 
#include <string>
#include <armadillo>
 
#include "equationofstate/equationofstatebase.hpp"
 
class BWRS : public EquationOfStateBase
{
public:
    BWRS(
        const arma::vec& composition = Composition::defaultComposition,
        const std::string& parameterSet = "Calsep"
    );
 
    static BWRS fromFilePaths(
        const arma::vec& composition = Composition::defaultComposition,
        const std::string& ABparameterFile = std::string(TRANSFLOW_RESOURCE_PATH) + "/equationofstate/bwrs/calsepABparameters.csv",
        const std::string& binaryInteractionTableFile = std::string(TRANSFLOW_RESOURCE_PATH) + "/equationofstate/bwrs/calsepBinaryInteraction.csv",
        const std::string& criticalProperties = "Starling"
    );
 
    virtual arma::vec evaluate(
            const double pressure,
            const double temperature) const override;
 
    virtual double calculateCompressibility(
            const double pressure,
            const double temperature) const override;
 
    virtual bool setComposition(const arma::vec& composition, const bool force = true) override;
 
    void loadCriticalProperties(const std::string name);
 
    void loadParametersAndCriticalProperties(const std::string parameterSet);
 
    void loadStarlingCriticalProperties();
 
    void loadJFHCriticalProperties();
 
    void loadCalsepCriticalProperties();
 
    void loadGasscoParameters();
 
    void loadCalsepParameters();
 
    void loadStarlingParameters();
 
    void loadParameterFiles(
            const std::string& ABparameterFile,
            const std::string& binaryInteractionTableFile);
 
    void enableConstantHeatCapacities();
 
    double findDensity(
            const double pressure,
            const double temperature,
            const double tolerance = 1e-4) const;
 
//    double getDerivativeOfZwrtDensityAtConstantTemperature(const double density, const double temperature) const;
 
    double getGasConstant() const { return m_R; } 
 
    double getMixtureCriticalPressure() const { return m_criticalPressureOfMixture; }
 
    double getMixtureCriticalTemperature() const { return m_criticalTemperatureOfMixture; }
 
protected:
 
    // parameters from Starling book
    //                      C1          C2          C3          iC4         nC4         iC5         nC5         C6          N2          CO2
 
    arma::vec m_Tc        = arma::vec({ 190.69,     305.39,     369.89,     408.13,     425.19,     460.37,     469.49,     507.29,     126.15,     304.15});           // critical temperature
    arma::vec m_rhoc      = arma::vec({ 1.00500e+4, 6.75659e+3, 4.99936e+3, 3.80118e+3, 3.92132e+3, 3.24694e+3, 3.21491e+3, 2.71673e+3, 1.10992e+4, 1.06379e+4});       // critical molar density [mol/m3]
    arma::vec m_w         = arma::vec({ 0.013,      0.1018,     0.157,      0.183,      0.197,      0.226,      0.252,      0.302,      0.035,      0.21});             // accentric factor
    arma::vec m_pc        = arma::vec({ 45.96,      48.839,     42.5,       36.48,      37.96,      33.81,      33.69,      27.34,      33.99,      73.825})*1e5;       // critical pressure
    arma::vec m_molarMass = arma::vec({ 16.042,     30.068,     44.094,     58.12,      58.12,      72.146,     72.146,     86.172,     28.016,     44.01});            // [g/mol]
    arma::vec m_expW = exp(-3.8*m_w);
 
    double m_R = 8.3160; // gas constant from Starling [m3 Pa / K mol]
 
    arma::mat m_binaryInteractionParameterTable;
 
    arma::vec m_Bi;
 
    arma::vec m_Ai;
 
    double m_A0; 
    double m_B0; 
    double m_C0; 
    double m_D0; 
    double m_E0; 
    double m_a; 
    double m_b; 
    double m_c; 
    double m_d; 
    double m_ALPHA; 
    double m_GAMMA; 
 
    double m_criticalPressureOfMixture; 
    double m_criticalTemperatureOfMixture; 
 
    bool m_useConstantHeatCapacities = false; 
 
    void calculateCoefficients();
 
    void findNonZeroComponents();
 
    arma::uvec m_indices; 
 
    double findMolarDensity(
            const double pressure,
            const double temperature,
            const double tolerance = 1e-4) const;
};