Method outerPQ

src/Backends/PCSAFT/PCSAFTBackend.cpp:2084–2343

}

double PCSAFTBackend::outerPQ(double t_guess, PCSAFTBackend& PCSAFT) {
    // Based on the algorithm proposed in H. A. J. Watson, M. Vikse, T. Gundersen, and P. I. Barton, “Reliable Flash Calculations: Part 1. Nonsmooth Inside-Out Algorithms,” Ind. Eng. Chem. Res., vol. 56, no. 4, pp. 960–973, Feb. 2017, doi: 10.1021/acs.iecr.6b03956.
    auto ncomp = N;  // number of components
    double TOL = 1e-8;
    int MAXITER = 200;

    // Define the residual to be driven to zero
    class SolverInnerResid : public FuncWrapper1D
    {
       public:
        PCSAFTBackend& PCSAFT;
        CoolPropDbl kb0;
        vector<CoolPropDbl> u;

        SolverInnerResid(PCSAFTBackend& PCSAFT, CoolPropDbl kb0, vector<CoolPropDbl> u) : PCSAFT(PCSAFT), kb0(kb0), u(std::move(u)) {}
        CoolPropDbl call(CoolPropDbl R) override {
            auto ncomp = PCSAFT.components.size();
            double error = 0;

            vector<double> pp(ncomp, 0);
            double L = 0;
            for (auto i = 0U; i < ncomp; i++) {
                if (!PCSAFT.ion_term || PCSAFT.components[i].getZ() == 0) {
                    pp[i] = PCSAFT.mole_fractions[i] / (1 - R + kb0 * R * exp(u[i]));
                    L += pp[i];
                } else {
                    L += PCSAFT.mole_fractions[i];
                }
            }
            L = (1 - R) * L;

            error = pow((L + PCSAFT._Q - 1), 2.);
            return error;
        };
    };

    double x_ions = 0.;  // overall mole fraction of ions in the system
    for (size_t i = 0; i < ncomp; i++) {
        if (PCSAFT.ion_term && PCSAFT.components[i].getZ() != 0) {
            x_ions += PCSAFT.mole_fractions[i];
        }
    }

    // initialize variables
    vector<double> k(ncomp, 0), u(ncomp, 0), kprime(ncomp, 0), uprime(ncomp, 0);
    double Tref = t_guess - 1;
    double Tprime = t_guess + 1;
    double t = t_guess;

    PCSAFT.SatL->_T = t;  // _T must be updated because the density calculation depends on it
    PCSAFT.SatV->_T = t;

    // calculate sigma for water, if it is present
    if (PCSAFT.water_present) {
        PCSAFT.components[water_idx].calc_water_sigma(t);
        PCSAFT.SatL->components[water_idx].calc_water_sigma(t);
        PCSAFT.SatV->components[water_idx].calc_water_sigma(t);
        PCSAFT.dielc =