compositionalflash.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_COMPOSITIONAL_FLASH_HH
#define DUMUX_COMPOSITIONAL_FLASH_HH
 
#include <dune/common/fvector.hh>
 
#include <dumux/material/fluidstates/pseudo1p2c.hh>
 
namespace Dumux {
 
template <class Scalar, class FluidSystem>
class CompositionalFlash
{
    using FluidState1p2c = PseudoOnePTwoCFluidState<Scalar, FluidSystem>;
 
    enum {  numPhases = FluidSystem::numPhases,
            numComponents = FluidSystem::numComponents
        };
 
    enum{
        phase0Idx = FluidSystem::phase0Idx,
        phase1Idx = FluidSystem::phase1Idx,
        comp0Idx = FluidSystem::comp0Idx,
        comp1Idx = FluidSystem::comp1Idx
    };
 
public:
    using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
    using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
 
 
    template<class FluidState>
    static void concentrationFlash2p2c(FluidState& fluidState,
                                       const Scalar Z0,
                                       const PhaseVector& phasePressure,
                                       const Scalar temperature)
    {
#ifndef NDEBUG
        // this solver can only handle fluid systems which
        // assume ideal mixtures of all fluids.
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            assert(FluidSystem::isIdealMixture(phaseIdx));
 
        }
#endif
 
        // set the temperature, pressure
        fluidState.setTemperature(temperature);
        fluidState.setPressure(phase0Idx, phasePressure[phase0Idx]);
        fluidState.setPressure(phase1Idx, phasePressure[phase1Idx]);
 
        // mole equilibrium ratios k for in case first phase is reference phase
        const Scalar k10 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx) * fluidState.pressure(phase0Idx)
                     / (FluidSystem::fugacityCoefficient(fluidState, phase1Idx, comp0Idx) * fluidState.pressure(phase1Idx));
        const Scalar k11 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx) * fluidState.pressure(phase0Idx)
                     / (FluidSystem::fugacityCoefficient(fluidState, phase1Idx, comp1Idx) * fluidState.pressure(phase1Idx));
 
        // get mole fraction from equilibrium constants
        fluidState.setMoleFraction(phase0Idx, comp0Idx, ((1. - k11) / (k10 - k11)));
        fluidState.setMoleFraction(phase1Idx, comp0Idx, (fluidState.moleFraction(phase0Idx,comp0Idx) * k10));
        fluidState.setMoleFraction(phase0Idx, comp1Idx, 1.0 - fluidState.moleFraction(phase0Idx,comp0Idx));
        fluidState.setMoleFraction(phase1Idx, comp1Idx, 1.0 - fluidState.moleFraction(phase1Idx,comp0Idx));
 
        // mass equilibrium ratios K for in case first phase is reference phase
        const Scalar K10 = fluidState.massFraction(phase1Idx, comp0Idx) / fluidState.massFraction(phase0Idx, comp0Idx);
        const Scalar K11 = (1. - fluidState.massFraction(phase1Idx, comp0Idx)) / (1. - fluidState.massFraction(phase0Idx, comp0Idx));
 
        // phase mass fraction Nu (ratio of phase mass to total phase mass) of first phase
        const Scalar Nu0 = 1. + ((Z0 * (K10 - 1.)) + ((1. - Z0) * (K11 - 1.))) / ((K11 - 1.) * (K10 - 1.));
 
        // an array of the phase mass fractions from which we will compute the saturations
        std::array<Scalar, 2> phaseMassFraction;
 
        // check phase presence
        if (Nu0 > 0. && Nu0 < 1.) // two phases present
            phaseMassFraction[phase0Idx] = Nu0;
        else if (Nu0 <= 0.) // only second phase present
        {
            phaseMassFraction[phase0Idx] = 0.0; // no first phase
            fluidState.setMassFraction(phase1Idx,comp0Idx, Z0); // assign complete mass dissolved into second phase
        }
        else // only first phase present
        {
            phaseMassFraction[phase0Idx] = 1.0; // no second phase
            fluidState.setMassFraction(phase0Idx, comp0Idx, Z0); // assign complete mass dissolved into first phase
        }
 
        // complete phase mass fractions
        phaseMassFraction[phase1Idx] = 1.0 - phaseMassFraction[phase0Idx];
 
        // get densities with correct composition
        fluidState.setDensity(phase0Idx, FluidSystem::density(fluidState, phase0Idx));
        fluidState.setDensity(phase1Idx, FluidSystem::density(fluidState, phase1Idx));
        fluidState.setMolarDensity(phase0Idx, FluidSystem::molarDensity(fluidState, phase0Idx));
        fluidState.setMolarDensity(phase1Idx, FluidSystem::molarDensity(fluidState, phase1Idx));
 
        fluidState.setViscosity(phase0Idx, FluidSystem::viscosity(fluidState, phase0Idx));
        fluidState.setViscosity(phase1Idx, FluidSystem::viscosity(fluidState, phase1Idx));
 
        Scalar sw = phaseMassFraction[phase0Idx] / fluidState.density(phase0Idx);
        sw /= (phaseMassFraction[phase0Idx] / fluidState.density(phase0Idx)
                    + phaseMassFraction[phase1Idx] / fluidState.density(phase1Idx));
        fluidState.setSaturation(phase0Idx, sw);
        fluidState.setSaturation(phase1Idx, 1.0-sw);
    }
 
    static void concentrationFlash1p2c(FluidState1p2c& fluidState, const Scalar& Z0,const Dune::FieldVector<Scalar,numPhases>
                                       phasePressure,const int presentPhaseIdx, const Scalar& temperature)
    {
        // set the temperature, pressure
        fluidState.setTemperature(temperature);
        fluidState.setPressure(phase0Idx, phasePressure[phase0Idx]);
        fluidState.setPressure(phase1Idx, phasePressure[phase1Idx]);
 
        fluidState.setPresentPhaseIdx(presentPhaseIdx);
        fluidState.setMassFraction(presentPhaseIdx,comp0Idx, Z0);
 
        // transform mass to mole fractions
        fluidState.setMoleFraction(presentPhaseIdx, comp0Idx, Z0 / FluidSystem::molarMass(comp0Idx)
                / (Z0 / FluidSystem::molarMass(comp0Idx) + (1. - Z0) / FluidSystem::molarMass(comp1Idx)));
 
        fluidState.setAverageMolarMass(presentPhaseIdx,
                fluidState.massFraction(presentPhaseIdx, comp0Idx) * FluidSystem::molarMass(comp0Idx)
                + fluidState.massFraction(presentPhaseIdx, comp1Idx) * FluidSystem::molarMass(comp1Idx));
 
        fluidState.setDensity(presentPhaseIdx, FluidSystem::density(fluidState, presentPhaseIdx));
        fluidState.setMolarDensity(presentPhaseIdx, FluidSystem::molarDensity(fluidState, presentPhaseIdx));
 
        fluidState.setViscosity(presentPhaseIdx, FluidSystem::viscosity(fluidState, presentPhaseIdx));
    }
 
 
    template<class FluidState>
    static void saturationFlash2p2c(FluidState& fluidState,
                                    const Scalar saturation,
                                    const PhaseVector& phasePressure,
                                    const Scalar temperature)
    {
#ifndef NDEBUG
        // this solver can only handle fluid systems which
        // assume ideal mixtures of all fluids.
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            assert(FluidSystem::isIdealMixture(phaseIdx));
 
        }
#endif
 
        // set the temperature, pressure
        fluidState.setTemperature(temperature);
        fluidState.setPressure(phase0Idx, phasePressure[phase0Idx]);
        fluidState.setPressure(phase1Idx, phasePressure[phase1Idx]);
 
        // mole equilibrium ratios k for in case first phase is reference phase
        const Scalar k10 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx) * fluidState.pressure(phase0Idx)
                     / (FluidSystem::fugacityCoefficient(fluidState, phase1Idx, comp0Idx) * fluidState.pressure(phase1Idx));
        const Scalar k11 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx) * fluidState.pressure(phase0Idx)
                     / (FluidSystem::fugacityCoefficient(fluidState, phase1Idx, comp1Idx) * fluidState.pressure(phase1Idx));
 
        // get mole fraction from equilibrium constants
        fluidState.setMoleFraction(phase0Idx,comp0Idx, ((1. - k11) / (k10 - k11)));
        fluidState.setMoleFraction(phase1Idx,comp0Idx, (fluidState.moleFraction(phase0Idx,comp0Idx) * k10));
        fluidState.setMoleFraction(phase0Idx, comp1Idx, 1.0 - fluidState.moleFraction(phase0Idx,comp0Idx));
        fluidState.setMoleFraction(phase1Idx, comp1Idx, 1.0 - fluidState.moleFraction(phase1Idx,comp0Idx));
 
        // get densities with correct composition
        fluidState.setDensity(phase0Idx, FluidSystem::density(fluidState, phase0Idx));
        fluidState.setDensity(phase1Idx, FluidSystem::density(fluidState, phase1Idx));
        fluidState.setMolarDensity(phase0Idx, FluidSystem::molarDensity(fluidState, phase0Idx));
        fluidState.setMolarDensity(phase1Idx, FluidSystem::molarDensity(fluidState, phase1Idx));
 
        fluidState.setViscosity(phase0Idx, FluidSystem::viscosity(fluidState, phase0Idx));
        fluidState.setViscosity(phase1Idx, FluidSystem::viscosity(fluidState, phase1Idx));
 
        // set saturation
        fluidState.setSaturation(phase0Idx, saturation);
        fluidState.setSaturation(phase1Idx, 1.0-saturation);
    }
};
 
} // end namespace Dumux
 
#endif