| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup ConstraintSolvers | ||
| 10 | * \brief Determines the pressures and saturations of all fluid phases | ||
| 11 | * given the total mass of all components. | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_COMPOSITIONAL_FLASH_HH | ||
| 14 | #define DUMUX_COMPOSITIONAL_FLASH_HH | ||
| 15 | |||
| 16 | #include <dune/common/fvector.hh> | ||
| 17 | |||
| 18 | #include <dumux/material/fluidstates/pseudo1p2c.hh> | ||
| 19 | |||
| 20 | namespace Dumux { | ||
| 21 | |||
| 22 | /*! | ||
| 23 | * \ingroup ConstraintSolvers | ||
| 24 | * \brief Flash calculation routines for compositional sequential models | ||
| 25 | * | ||
| 26 | * Routines for isothermal and isobaric 2p2c and 1p2c flash, assuming an ideal mixture. | ||
| 27 | * The flash assumes that the fugacities of a component \f$ \kappa \f$ in each phase | ||
| 28 | * are the same. The fugacity is defined as: | ||
| 29 | * | ||
| 30 | * \f$ f^\kappa = \Phi^\kappa_\alpha(T_\alpha, p_\alpha) p_\alpha x^\kappa_\alpha\; \f$, | ||
| 31 | * | ||
| 32 | * where \f$ Phi^\kappa_\alpha \f$ is a fixed fugacity coefficient independent of the phase's composition | ||
| 33 | * (ideal mixture), \f$ T_\alpha \f$ is a fixed temperature, and \f$ p_\alpha \f$ a fixed phase pressure. | ||
| 34 | * From the equality of fugacities, the mole (and mass) fractions \f$ x^\kappa_\alpha \f$ | ||
| 35 | * in equilibrium are calculated. | ||
| 36 | */ | ||
| 37 | template <class Scalar, class FluidSystem> | ||
| 38 | class CompositionalFlash | ||
| 39 | { | ||
| 40 | using FluidState1p2c = PseudoOnePTwoCFluidState<Scalar, FluidSystem>; | ||
| 41 | |||
| 42 | enum { numPhases = FluidSystem::numPhases, | ||
| 43 | numComponents = FluidSystem::numComponents | ||
| 44 | }; | ||
| 45 | |||
| 46 | enum{ | ||
| 47 | phase0Idx = FluidSystem::phase0Idx, | ||
| 48 | phase1Idx = FluidSystem::phase1Idx, | ||
| 49 | comp0Idx = FluidSystem::comp0Idx, | ||
| 50 | comp1Idx = FluidSystem::comp1Idx | ||
| 51 | }; | ||
| 52 | |||
| 53 | public: | ||
| 54 | using ComponentVector = Dune::FieldVector<Scalar, numComponents>; | ||
| 55 | using PhaseVector = Dune::FieldVector<Scalar, numPhases>; | ||
| 56 | |||
| 57 | |||
| 58 | /*! | ||
| 59 | * \name Concentration flash for a given feed fraction | ||
| 60 | * \brief 2p2c Flash for constant p & T if concentration (feed mass fraction) is given. | ||
| 61 | * | ||
| 62 | * This flash uses the Rachford-Rice equation: | ||
| 63 | * Rachford Jr, H. H., & Rice, J. D. (1952). | ||
| 64 | * Procedure for use of electronic digital computers in calculating flash vaporization | ||
| 65 | * hydrocarbon equilibrium. Journal of Petroleum Technology, 4(10), 19-3. | ||
| 66 | * | ||
| 67 | * Routine goes as follows: | ||
| 68 | * - determination of the equilibrium constants from the fluid system | ||
| 69 | * - determination of maximum solubilities (mole fractions) according to phase pressures | ||
| 70 | * - comparison with phase mass fraction Nu (from Rachford-Rice equation) | ||
| 71 | * to determine phase presence => actual mole and mass fractions | ||
| 72 | * - complete fluid state | ||
| 73 | * \param fluidState The sequential fluid State | ||
| 74 | * \param Z0 Feed mass fraction: Mass of first component per total mass \f$\mathrm{[-]}\f$ | ||
| 75 | * \param phasePressure Vector holding the pressure \f$\mathrm{[Pa]}\f$ | ||
| 76 | * \param temperature Temperature \f$\mathrm{[K]}\f$ | ||
| 77 | */ | ||
| 78 | template<class FluidState> | ||
| 79 | 7 | static void concentrationFlash2p2c(FluidState& fluidState, | |
| 80 | const Scalar Z0, | ||
| 81 | const PhaseVector& phasePressure, | ||
| 82 | const Scalar temperature) | ||
| 83 | { | ||
| 84 | #ifndef NDEBUG | ||
| 85 | // this solver can only handle fluid systems which | ||
| 86 | // assume ideal mixtures of all fluids. | ||
| 87 | 7 | for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { | |
| 88 | assert(FluidSystem::isIdealMixture(phaseIdx)); | ||
| 89 | |||
| 90 | } | ||
| 91 | #endif | ||
| 92 | |||
| 93 | // set the temperature, pressure | ||
| 94 | 7 | fluidState.setTemperature(temperature); | |
| 95 | 14 | fluidState.setPressure(phase0Idx, phasePressure[phase0Idx]); | |
| 96 | 14 | fluidState.setPressure(phase1Idx, phasePressure[phase1Idx]); | |
| 97 | |||
| 98 | // mole equilibrium ratios k for in case first phase is reference phase | ||
| 99 | 7 | const Scalar k10 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx) * fluidState.pressure(phase0Idx) | |
| 100 | 7 | / (FluidSystem::fugacityCoefficient(fluidState, phase1Idx, comp0Idx) * fluidState.pressure(phase1Idx)); | |
| 101 | 7 | const Scalar k11 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx) * fluidState.pressure(phase0Idx) | |
| 102 | 7 | / (FluidSystem::fugacityCoefficient(fluidState, phase1Idx, comp1Idx) * fluidState.pressure(phase1Idx)); | |
| 103 | |||
| 104 | // get mole fraction from equilibrium constants | ||
| 105 | 7 | fluidState.setMoleFraction(phase0Idx, comp0Idx, ((1. - k11) / (k10 - k11))); | |
| 106 | 14 | fluidState.setMoleFraction(phase1Idx, comp0Idx, (fluidState.moleFraction(phase0Idx,comp0Idx) * k10)); | |
| 107 | 14 | fluidState.setMoleFraction(phase0Idx, comp1Idx, 1.0 - fluidState.moleFraction(phase0Idx,comp0Idx)); | |
| 108 | 14 | fluidState.setMoleFraction(phase1Idx, comp1Idx, 1.0 - fluidState.moleFraction(phase1Idx,comp0Idx)); | |
| 109 | |||
| 110 | // mass equilibrium ratios K for in case first phase is reference phase | ||
| 111 |
4/4✓ Branch 0 taken 5 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 2 times.
|
14 | const Scalar K10 = fluidState.massFraction(phase1Idx, comp0Idx) / fluidState.massFraction(phase0Idx, comp0Idx); |
| 112 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 2 times.
|
7 | const Scalar K11 = (1. - fluidState.massFraction(phase1Idx, comp0Idx)) / (1. - fluidState.massFraction(phase0Idx, comp0Idx)); |
| 113 | |||
| 114 | // phase mass fraction Nu (ratio of phase mass to total phase mass) of first phase | ||
| 115 | 7 | const Scalar Nu0 = 1. + ((Z0 * (K10 - 1.)) + ((1. - Z0) * (K11 - 1.))) / ((K11 - 1.) * (K10 - 1.)); | |
| 116 | |||
| 117 | // an array of the phase mass fractions from which we will compute the saturations | ||
| 118 | std::array<Scalar, 2> phaseMassFraction; | ||
| 119 | |||
| 120 | // check phase presence | ||
| 121 |
4/4✓ Branch 0 taken 5 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 3 times.
|
7 | if (Nu0 > 0. && Nu0 < 1.) // two phases present |
| 122 | phaseMassFraction[phase0Idx] = Nu0; | ||
| 123 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
|
4 | else if (Nu0 <= 0.) // only second phase present |
| 124 | { | ||
| 125 | 2 | phaseMassFraction[phase0Idx] = 0.0; // no first phase | |
| 126 | 2 | fluidState.setMassFraction(phase1Idx,comp0Idx, Z0); // assign complete mass dissolved into second phase | |
| 127 | } | ||
| 128 | else // only first phase present | ||
| 129 | { | ||
| 130 | 2 | phaseMassFraction[phase0Idx] = 1.0; // no second phase | |
| 131 | 2 | fluidState.setMassFraction(phase0Idx, comp0Idx, Z0); // assign complete mass dissolved into first phase | |
| 132 | } | ||
| 133 | |||
| 134 | // complete phase mass fractions | ||
| 135 | 14 | phaseMassFraction[phase1Idx] = 1.0 - phaseMassFraction[phase0Idx]; | |
| 136 | |||
| 137 | // get densities with correct composition | ||
| 138 | 7 | fluidState.setDensity(phase0Idx, FluidSystem::density(fluidState, phase0Idx)); | |
| 139 | 7 | fluidState.setDensity(phase1Idx, FluidSystem::density(fluidState, phase1Idx)); | |
| 140 | 7 | fluidState.setMolarDensity(phase0Idx, FluidSystem::molarDensity(fluidState, phase0Idx)); | |
| 141 | 14 | fluidState.setMolarDensity(phase1Idx, FluidSystem::molarDensity(fluidState, phase1Idx)); | |
| 142 | |||
| 143 | 7 | fluidState.setViscosity(phase0Idx, FluidSystem::viscosity(fluidState, phase0Idx)); | |
| 144 | 7 | fluidState.setViscosity(phase1Idx, FluidSystem::viscosity(fluidState, phase1Idx)); | |
| 145 | |||
| 146 | 21 | Scalar sw = phaseMassFraction[phase0Idx] / fluidState.density(phase0Idx); | |
| 147 | 28 | sw /= (phaseMassFraction[phase0Idx] / fluidState.density(phase0Idx) | |
| 148 | 21 | + phaseMassFraction[phase1Idx] / fluidState.density(phase1Idx)); | |
| 149 | 14 | fluidState.setSaturation(phase0Idx, sw); | |
| 150 | 14 | fluidState.setSaturation(phase1Idx, 1.0-sw); | |
| 151 | 7 | } | |
| 152 | |||
| 153 | /*! | ||
| 154 | * \brief The simplest possible update routine for 1p2c "flash" calculations | ||
| 155 | * | ||
| 156 | * Routine goes as follows: | ||
| 157 | * - check if we are in single phase condition | ||
| 158 | * - assign total concentration to the present phase | ||
| 159 | * - complete fluid state | ||
| 160 | * | ||
| 161 | * \param fluidState The sequential fluid state | ||
| 162 | * \param Z0 Feed mass fraction: Mass of first component per total mass \f$\mathrm{[-]}\f$ | ||
| 163 | * \param phasePressure Vector holding the pressure \f$\mathrm{[Pa]}\f$ | ||
| 164 | * \param presentPhaseIdx Subdomain Index = Indication which phase is present | ||
| 165 | * \param temperature Temperature \f$\mathrm{[K]}\f$ | ||
| 166 | */ | ||
| 167 | 2 | static void concentrationFlash1p2c(FluidState1p2c& fluidState, const Scalar& Z0,const Dune::FieldVector<Scalar,numPhases> | |
| 168 | phasePressure,const int presentPhaseIdx, const Scalar& temperature) | ||
| 169 | { | ||
| 170 | // set the temperature, pressure | ||
| 171 | 2 | fluidState.setTemperature(temperature); | |
| 172 | 4 | fluidState.setPressure(phase0Idx, phasePressure[phase0Idx]); | |
| 173 | 4 | fluidState.setPressure(phase1Idx, phasePressure[phase1Idx]); | |
| 174 | |||
| 175 | 2 | fluidState.setPresentPhaseIdx(presentPhaseIdx); | |
| 176 | 2 | fluidState.setMassFraction(presentPhaseIdx,comp0Idx, Z0); | |
| 177 | |||
| 178 | // transform mass to mole fractions | ||
| 179 | 6 | fluidState.setMoleFraction(presentPhaseIdx, comp0Idx, Z0 / FluidSystem::molarMass(comp0Idx) | |
| 180 | 2 | / (Z0 / FluidSystem::molarMass(comp0Idx) + (1. - Z0) / FluidSystem::molarMass(comp1Idx))); | |
| 181 | |||
| 182 | 4 | fluidState.setAverageMolarMass(presentPhaseIdx, | |
| 183 | 2 | fluidState.massFraction(presentPhaseIdx, comp0Idx) * FluidSystem::molarMass(comp0Idx) | |
| 184 | 4 | + fluidState.massFraction(presentPhaseIdx, comp1Idx) * FluidSystem::molarMass(comp1Idx)); | |
| 185 | |||
| 186 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
|
2 | fluidState.setDensity(presentPhaseIdx, FluidSystem::density(fluidState, presentPhaseIdx)); |
| 187 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
|
2 | fluidState.setMolarDensity(presentPhaseIdx, FluidSystem::molarDensity(fluidState, presentPhaseIdx)); |
| 188 | |||
| 189 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
|
2 | fluidState.setViscosity(presentPhaseIdx, FluidSystem::viscosity(fluidState, presentPhaseIdx)); |
| 190 | 2 | } | |
| 191 | //@} | ||
| 192 | |||
| 193 | //@{ | ||
| 194 | /*! \name Saturation flash for a given saturation (e.g. at boundary) | ||
| 195 | * \brief 2p2c flash for constant p & T if the saturation instead of concentration (feed mass fraction) is known. | ||
| 196 | * | ||
| 197 | * Routine goes as follows: | ||
| 198 | * - determination of the equilibrium constants from the fluid system | ||
| 199 | * - determination of maximum solubilities (mole fractions) according to phase pressures | ||
| 200 | * - complete fluid state | ||
| 201 | * \param fluidState The sequential fluid state | ||
| 202 | * \param saturation Saturation of phase 1 \f$\mathrm{[-]}\f$ | ||
| 203 | * \param phasePressure Vector holding the pressure \f$\mathrm{[Pa]}\f$ | ||
| 204 | * \param temperature Temperature \f$\mathrm{[K]}\f$ | ||
| 205 | */ | ||
| 206 | template<class FluidState> | ||
| 207 | 1 | static void saturationFlash2p2c(FluidState& fluidState, | |
| 208 | const Scalar saturation, | ||
| 209 | const PhaseVector& phasePressure, | ||
| 210 | const Scalar temperature) | ||
| 211 | { | ||
| 212 | #ifndef NDEBUG | ||
| 213 | // this solver can only handle fluid systems which | ||
| 214 | // assume ideal mixtures of all fluids. | ||
| 215 | 1 | for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { | |
| 216 | assert(FluidSystem::isIdealMixture(phaseIdx)); | ||
| 217 | |||
| 218 | } | ||
| 219 | #endif | ||
| 220 | |||
| 221 | // set the temperature, pressure | ||
| 222 | 1 | fluidState.setTemperature(temperature); | |
| 223 | 2 | fluidState.setPressure(phase0Idx, phasePressure[phase0Idx]); | |
| 224 | 2 | fluidState.setPressure(phase1Idx, phasePressure[phase1Idx]); | |
| 225 | |||
| 226 | // mole equilibrium ratios k for in case first phase is reference phase | ||
| 227 | 1 | const Scalar k10 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx) * fluidState.pressure(phase0Idx) | |
| 228 | 1 | / (FluidSystem::fugacityCoefficient(fluidState, phase1Idx, comp0Idx) * fluidState.pressure(phase1Idx)); | |
| 229 | 1 | const Scalar k11 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx) * fluidState.pressure(phase0Idx) | |
| 230 | 1 | / (FluidSystem::fugacityCoefficient(fluidState, phase1Idx, comp1Idx) * fluidState.pressure(phase1Idx)); | |
| 231 | |||
| 232 | // get mole fraction from equilibrium constants | ||
| 233 | 1 | fluidState.setMoleFraction(phase0Idx,comp0Idx, ((1. - k11) / (k10 - k11))); | |
| 234 | 2 | fluidState.setMoleFraction(phase1Idx,comp0Idx, (fluidState.moleFraction(phase0Idx,comp0Idx) * k10)); | |
| 235 | 2 | fluidState.setMoleFraction(phase0Idx, comp1Idx, 1.0 - fluidState.moleFraction(phase0Idx,comp0Idx)); | |
| 236 | 2 | fluidState.setMoleFraction(phase1Idx, comp1Idx, 1.0 - fluidState.moleFraction(phase1Idx,comp0Idx)); | |
| 237 | |||
| 238 | // get densities with correct composition | ||
| 239 | 1 | fluidState.setDensity(phase0Idx, FluidSystem::density(fluidState, phase0Idx)); | |
| 240 | 1 | fluidState.setDensity(phase1Idx, FluidSystem::density(fluidState, phase1Idx)); | |
| 241 | 1 | fluidState.setMolarDensity(phase0Idx, FluidSystem::molarDensity(fluidState, phase0Idx)); | |
| 242 | 2 | fluidState.setMolarDensity(phase1Idx, FluidSystem::molarDensity(fluidState, phase1Idx)); | |
| 243 | |||
| 244 | 1 | fluidState.setViscosity(phase0Idx, FluidSystem::viscosity(fluidState, phase0Idx)); | |
| 245 | 1 | fluidState.setViscosity(phase1Idx, FluidSystem::viscosity(fluidState, phase1Idx)); | |
| 246 | |||
| 247 | // set saturation | ||
| 248 | 2 | fluidState.setSaturation(phase0Idx, saturation); | |
| 249 | 2 | fluidState.setSaturation(phase1Idx, 1.0-saturation); | |
| 250 | 1 | } | |
| 251 | //@} | ||
| 252 | }; | ||
| 253 | |||
| 254 | } // end namespace Dumux | ||
| 255 | |||
| 256 | #endif | ||
| 257 |