UNIFAC.cpp

Go to the documentation of this file.

#include "UNIFAC.h"
 
void UNIFAC::UNIFACMixture::set_interaction_parameters() {
    for (const auto& isgi : unique_groups) {
        for (const auto& jsgi : unique_groups) {
            if (jsgi >= isgi) {
                continue;
            }
            int mgi1 = static_cast<int>(m_sgi_to_mgi.find(isgi)->second);
            int mgi2 = static_cast<int>(m_sgi_to_mgi.find(jsgi)->second);
            // Insert in normal order
            std::pair<std::pair<int, int>, UNIFACLibrary::InteractionParameters> m_pair(std::pair<int, int>(mgi1, mgi2),
                                                                                        library.get_interaction_parameters(mgi1, mgi2));
            interaction.insert(m_pair);
            // Insert in backwards order
            if (mgi1 != mgi2) {
                std::pair<std::pair<int, int>, UNIFACLibrary::InteractionParameters> m_pair(std::pair<int, int>(mgi2, mgi1),
                                                                                            library.get_interaction_parameters(mgi2, mgi1));
                interaction.insert(m_pair);
            }
        }
    }
}
 
void UNIFAC::UNIFACMixture::set_interaction_parameter(const std::size_t mgi1, const std::size_t mgi2, const std::string& parameter,
                                                      const double value) {
    if (parameter == "aij") {
        this->interaction[std::pair<int, int>(static_cast<int>(mgi1), static_cast<int>(mgi2))].a_ij = value;
    } else if (parameter == "bij") {
        this->interaction[std::pair<int, int>(static_cast<int>(mgi1), static_cast<int>(mgi2))].b_ij = value;
    } else if (parameter == "cij") {
        this->interaction[std::pair<int, int>(static_cast<int>(mgi1), static_cast<int>(mgi2))].c_ij = value;
    } else {
        throw CoolProp::ValueError(format("I don't know what to do with parameter [%s]", parameter.c_str()));
    }
}
double UNIFAC::UNIFACMixture::get_interaction_parameter(const std::size_t mgi1, const std::size_t mgi2, const std::string& parameter) {
    auto it = this->interaction.find(std::pair<int, int>(static_cast<int>(mgi1), static_cast<int>(mgi2)));
    if (it == this->interaction.end()) {
        throw CoolProp::ValueError(format("Unable to match mgi-mgi pair: [%d,%d]", static_cast<int>(mgi1), static_cast<int>(mgi1)));
    } else {
        if (parameter == "aij") {
            return it->second.a_ij;
        } else if (parameter == "bij") {
            return it->second.b_ij;
        } else if (parameter == "cij") {
            return it->second.c_ij;
        } else {
            throw CoolProp::ValueError(format("I don't know what to do with parameter [%s]", parameter.c_str()));
        }
    }
}
 
void UNIFAC::UNIFACMixture::set_mole_fractions(const std::vector<double>& z) {
    //    // If the vector fractions are the same as last ones, don't do anything and return
    //    if (!mole_fractions.empty() && maxvectordiff(z, mole_fractions) < 1e-15){
    //        return;
    //    }
    this->mole_fractions = z;
    if (this->N != z.size()) {
        throw CoolProp::ValueError("Size of molar fraction do not match number of components.");
    }
 
    std::map<std::size_t, double>&Xg = m_Xg, &thetag = m_thetag;
    Xg.clear();
    thetag.clear();
 
    // Iterate over the fluids
    double X_summer = 0;
    for (std::size_t i = 0; i < this->mole_fractions.size(); ++i) {
        X_summer += this->mole_fractions[i] * pure_data[i].group_count;
    }
    for (const auto& sgi : unique_groups) {
        double X = 0;
        // Iterate over the fluids
        for (std::size_t i = 0; i < this->mole_fractions.size(); ++i) {
            X += this->mole_fractions[i] * group_count(i, sgi);
        }
        Xg.emplace(sgi, X);
    }
    for (auto& [key, val] : Xg) {
        val /= X_summer;
        //printf("X_{%d}: %g\n", key, val);
    }
    double theta_summer = 0;
    for (const auto& sgi : unique_groups) {
        double cont = Xg.find(sgi)->second * m_Q.find(sgi)->second;
        theta_summer += cont;
        thetag.emplace(sgi, cont);
    }
    for (auto& [key, val] : thetag) {
        val /= theta_summer;
        //printf("theta_{%d}: %g\n", key, val);
    }
}
 
double UNIFAC::UNIFACMixture::Psi(std::size_t sgi1, std::size_t sgi2) const {
 
    if (this->interaction.size() == 0) {
        throw CoolProp::ValueError("interaction parameters for UNIFAC not yet set");
    }
    std::size_t mgi1 = m_sgi_to_mgi.find(sgi1)->second;
    std::size_t mgi2 = m_sgi_to_mgi.find(sgi2)->second;
    if (mgi1 == mgi2) {
        return 1;
    } else {
        auto it = this->interaction.find(std::pair<int, int>(static_cast<int>(mgi1), static_cast<int>(mgi2)));
        if (it != this->interaction.end()) {
            return exp(-(it->second.a_ij / this->m_T + it->second.b_ij + it->second.c_ij * this->m_T));
        } else {
            throw CoolProp::ValueError(
              format("Could not match mgi[%d]-mgi[%d] interaction in UNIFAC", static_cast<int>(mgi1), static_cast<int>(mgi2)));
        }
    }
}
 
std::size_t UNIFAC::UNIFACMixture::group_count(std::size_t i, std::size_t sgi) const {
    const UNIFACLibrary::Component& c = components[i];
    for (const auto& cg : c.groups) {
        if (cg.group.sgi == sgi) {
            return cg.count;
        }
    }
    return 0;
}
 
double UNIFAC::UNIFACMixture::theta_pure(std::size_t i, std::size_t sgi) const {
    return pure_data[i].theta.find(sgi)->second;
}
 
void UNIFAC::UNIFACMixture::set_temperature(const double T) {
    //    // Check whether you are using exactly the same temperature as last time
    //    if (static_cast<bool>(_T) && std::abs(static_cast<double>(_T) - T) < 1e-15 {
    //        //
    //        return;
    //    }
    this->m_T = T;
 
    if (this->mole_fractions.empty()) {
        throw CoolProp::ValueError("mole fractions must be set before calling set_temperature");
    }
 
    // Compute Psi once for the different calls
    for (const auto& k : unique_groups) {
        for (const auto& m : unique_groups) {
            Psi_[std::pair<std::size_t, std::size_t>(k, m)] = Psi(k, m);
        }
    }
 
    for (std::size_t i = 0; i < this->mole_fractions.size(); ++i) {
        const UNIFACLibrary::Component& c = components[i];
        for (std::size_t k = 0; k < c.groups.size(); ++k) {
            double Q = c.groups[k].group.Q_k;
            int sgik = c.groups[k].group.sgi;
            double sum1 = 0;
            for (const auto& group : c.groups) {
                int sgim = group.group.sgi;
                sum1 += theta_pure(i, sgim) * Psi_.find(std::pair<std::size_t, std::size_t>(sgim, sgik))->second;
            }
            // sum1 > 0: c.groups is the set of groups present in component i, theta_pure
            // returns the (positive) surface-area fraction of each, and Psi = exp(-a/T) > 0.
            // cppcheck cannot prove this symbolically.
            // cppcheck-suppress invalidFunctionArg
            double s = 1 - log(sum1);
            for (std::size_t m = 0; m < c.groups.size(); ++m) {
                int sgim = c.groups[m].group.sgi;
                double sum2 = 0;
                for (const auto& group : c.groups) {
                    int sgin = group.group.sgi;
                    sum2 += theta_pure(i, sgin) * Psi_.find(std::pair<std::size_t, std::size_t>(sgin, sgim))->second;
                }
                s -= theta_pure(i, sgim) * Psi_.find(std::pair<std::size_t, std::size_t>(sgik, sgim))->second / sum2;
            }
            pure_data[i].lnGamma[sgik] = Q * s;
            //printf("ln(Gamma)^(%d)_{%d}: %g\n", static_cast<int>(i + 1), sgik, Q*s);
        }
    }
 
    std::map<std::size_t, double>&thetag = m_thetag, &lnGammag = m_lnGammag;
    lnGammag.clear();
 
    for (const auto& ksgi : unique_groups) {
        double sum1 = 0;
        for (const auto& msgi : unique_groups) {
            sum1 += thetag.find(msgi)->second * Psi_.find(std::pair<std::size_t, std::size_t>(msgi, ksgi))->second;
        }
        double s = 1 - log(sum1);
        for (const auto& msgi : unique_groups) {
            double sum3 = 0;
            for (const auto& nsgi : unique_groups) {
                sum3 += thetag.find(nsgi)->second * Psi_.find(std::pair<std::size_t, std::size_t>(nsgi, msgi))->second;
            }
            s -= thetag.find(msgi)->second * Psi_.find(std::pair<std::size_t, std::size_t>(ksgi, msgi))->second / sum3;
        }
        lnGammag.emplace(ksgi, m_Q.find(ksgi)->second * s);
        //printf("log(Gamma)_{%d}: %g\n", itk->sgi, itk->Q_k*s);
    }
    _T = m_T;
}
double UNIFAC::UNIFACMixture::ln_gamma_R(const double tau, std::size_t i, std::size_t itau) {
    if (itau == 0) {
        set_temperature(T_r / tau);
        double summer = 0;
        for (const auto& sgi : unique_groups) {
            std::size_t count = group_count(i, sgi);
            if (count > 0) {
                summer += count * (m_lnGammag.find(sgi)->second - pure_data[i].lnGamma.find(sgi)->second);
            }
        }
        //printf("log(gamma)_{%d}: %g\n", i+1, summer);
        return summer;
    } else {
        double dtau = 0.01 * tau;
        return (ln_gamma_R(tau + dtau, i, itau - 1) - ln_gamma_R(tau - dtau, i, itau - 1)) / (2 * dtau);
    }
}
void UNIFAC::UNIFACMixture::activity_coefficients(double tau, const std::vector<double>& z, std::vector<double>& gamma) {
    if (this->N != z.size()) {
        throw CoolProp::ValueError("Size of molar fraction do not match number of components.");
    }
    std::vector<double> r(N), q(N), l(N), phi(N), theta(N), ln_Gamma_C(N);
    double summerzr = 0, summerzq = 0, summerzl = 0;
    for (std::size_t i = 0; i < N; ++i) {
        double summerr = 0, summerq = 0;
        const UNIFACLibrary::Component& c = components[i];
        for (const auto& cg : c.groups) {
            summerr += cg.count * cg.group.R_k;
            summerq += cg.count * cg.group.Q_k;
        }
        r[i] = summerr;
        q[i] = summerq;
        summerzr += z[i] * r[i];
        summerzq += z[i] * q[i];
    }
    for (std::size_t i = 0; i < N; ++i) {
        phi[i] = z[i] * r[i] / summerzr;
        theta[i] = z[i] * q[i] / summerzq;
        l[i] = 10.0 / 2.0 * (r[i] - q[i]) - (r[i] - 1);
        summerzl += z[i] * l[i];
    }
    for (std::size_t i = 0; i < N; ++i) {
        ln_Gamma_C[i] = log(phi[i] / z[i]) + 10.0 / 2.0 * q[i] * log(theta[i] / phi[i]) + l[i] - phi[i] / z[i] * summerzl;
        gamma[i] = exp(ln_gamma_R(tau, i, 0) + ln_Gamma_C[i]);
    }
}
 
void UNIFAC::UNIFACMixture::add_component(const UNIFACLibrary::Component& comp) {
    components.push_back(comp);
    for (const auto& cg : comp.groups) {
        m_sgi_to_mgi.emplace(cg.group.sgi, cg.group.mgi);
    }
}
 
void UNIFAC::UNIFACMixture::set_components(const std::string& identifier_type, const std::vector<std::string>& identifiers) {
    components.clear();
    N = identifiers.size();
    if (identifier_type == "name") {
        // Iterate over the provided names
        for (const auto& id : identifiers) {
            // Get and add the component
            UNIFACLibrary::Component c = library.get_component("name", id);
            add_component(c);
        }
    } else {
        throw CoolProp::ValueError("Cannot understand identifier_type");
    }
    set_pure_data();
}
 
void UNIFAC::UNIFACMixture::set_pure_data() {
    pure_data.clear();
    unique_groups.clear();
    m_Q.clear();
    for (std::size_t i = 0; i < N; ++i) {
        const UNIFACLibrary::Component& c = components[i];
        ComponentData cd;
        double summerxq = 0;
        cd.group_count = 0;
        for (const auto& cg : c.groups) {
            auto x = static_cast<double>(cg.count);
            auto theta = static_cast<double>(cg.count * cg.group.Q_k);
            cd.X.emplace(cg.group.sgi, x);
            cd.theta.emplace(cg.group.sgi, theta);
            cd.group_count += cg.count;
            summerxq += x * cg.group.Q_k;
            unique_groups.insert(cg.group.sgi);
            m_Q.emplace(cg.group.sgi, cg.group.Q_k);
        }
        for (auto& [key, val] : cd.X) {
            val /= cd.group_count;
            //printf("X^(%d)_{%d}: %g\n", static_cast<int>(i + 1), static_cast<int>(key), val);
        }
        for (auto& [key, val] : cd.theta) {
            val /= summerxq;
            //printf("theta^(%d)_{%d}: %g\n", static_cast<int>(i+1), static_cast<int>(key), val);
        }
        pure_data.push_back(cd);
    }
}
 
void UNIFAC::UNIFACMixture::set_Q_k(const size_t sgi, const double value) {
    for (std::size_t i = 0; i < N; ++i) {
        for (auto& group : components[i].groups) {
            if (group.group.sgi == sgi) {
                group.group.Q_k = value;
            }
        }
    }
 
    set_pure_data();
}
 
double UNIFAC::UNIFACMixture::get_Q_k(const size_t sgi) const {
    for (std::size_t i = 0; i < N; ++i) {
        for (const auto& group : components[i].groups) {
            if (group.group.sgi == sgi) {
                return group.group.Q_k;
            }
        }
    }
    throw CoolProp::ValueError("Could not get Q_k");
}