| 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 FluidSystems | ||
| 10 | * \brief @copybrief Dumux::FluidSystems::BrineCO2 | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_BRINE_CO2_FLUID_SYSTEM_HH | ||
| 13 | #define DUMUX_BRINE_CO2_FLUID_SYSTEM_HH | ||
| 14 | |||
| 15 | #include <type_traits> | ||
| 16 | |||
| 17 | #include <dune/common/exceptions.hh> | ||
| 18 | |||
| 19 | #include <dumux/common/parameters.hh> | ||
| 20 | #include <dumux/material/idealgas.hh> | ||
| 21 | #include <dumux/material/fluidsystems/base.hh> | ||
| 22 | #include <dumux/material/fluidsystems/brine.hh> | ||
| 23 | #include <dumux/material/fluidstates/adapter.hh> | ||
| 24 | |||
| 25 | #include <dumux/material/components/brine.hh> | ||
| 26 | #include <dumux/material/components/co2.hh> | ||
| 27 | #include <dumux/material/components/tabulatedcomponent.hh> | ||
| 28 | |||
| 29 | #include <dumux/material/binarycoefficients/brine_co2.hh> | ||
| 30 | |||
| 31 | #include <dumux/io/name.hh> | ||
| 32 | |||
| 33 | namespace Dumux::FluidSystems::Detail { | ||
| 34 | |||
| 35 | /*! | ||
| 36 | * \brief Class that exports some indices that should | ||
| 37 | * be provided by the BrineCO2 fluid system. | ||
| 38 | * The indices are chosen dependent on the policy, | ||
| 39 | * i.e. if a simplified pseudo component Brine is | ||
| 40 | * used or salt is considered an individual component. | ||
| 41 | */ | ||
| 42 | template<bool useConstantSalinity> | ||
| 43 | struct BrineCO2Indices; | ||
| 44 | |||
| 45 | /*! | ||
| 46 | * \brief Specialization for the case of brine being | ||
| 47 | * a pseudo component with a constant salinity. | ||
| 48 | * \note This specialization exports brine as component | ||
| 49 | */ | ||
| 50 | template<> | ||
| 51 | struct BrineCO2Indices<true> | ||
| 52 | { | ||
| 53 | static constexpr int BrineIdx = 0; | ||
| 54 | }; | ||
| 55 | |||
| 56 | /*! | ||
| 57 | * \brief Specialization for the case of brine being | ||
| 58 | * a fluid system with NaCl as individual component. | ||
| 59 | * \note This specialization exports both water and NaCl as components | ||
| 60 | */ | ||
| 61 | template<> | ||
| 62 | struct BrineCO2Indices<false> | ||
| 63 | { | ||
| 64 | static constexpr int H2OIdx = 0; | ||
| 65 | static constexpr int NaClIdx = 2; | ||
| 66 | static constexpr int comp2Idx = 2; | ||
| 67 | }; | ||
| 68 | |||
| 69 | } // end namespace Dumux::FluidSystems::Detail | ||
| 70 | |||
| 71 | |||
| 72 | namespace Dumux::FluidSystems { | ||
| 73 | |||
| 74 | /*! | ||
| 75 | * \ingroup FluidSystems | ||
| 76 | * \brief Default policy for the Brine-CO2 fluid system | ||
| 77 | */ | ||
| 78 | template<bool salinityIsConstant, bool fastButSimplifiedRelations = false> | ||
| 79 | struct BrineCO2DefaultPolicy | ||
| 80 | { | ||
| 81 | static constexpr bool useConstantSalinity() { return salinityIsConstant; } | ||
| 82 | static constexpr bool useCO2GasDensityAsGasMixtureDensity() { return fastButSimplifiedRelations; } | ||
| 83 | }; | ||
| 84 | |||
| 85 | /*! | ||
| 86 | * \ingroup FluidSystems | ||
| 87 | * \brief A compositional fluid with brine (H2O & NaCl) and carbon dioxide as | ||
| 88 | * components in both the liquid and the gas (supercritical) phase. | ||
| 89 | * | ||
| 90 | * \note Depending on the chosen policy, the salinity is assumed to be constant | ||
| 91 | * (in which case Brine is used as a pseudo component) or salt (here NaCl) | ||
| 92 | * is considered as an individual component. | ||
| 93 | * \note This implementation always assumes NaCl stays in the liquid phase. | ||
| 94 | */ | ||
| 95 | template< class Scalar, | ||
| 96 | class CO2Component, | ||
| 97 | class H2OType = Components::TabulatedComponent<Components::H2O<Scalar>>, | ||
| 98 | class Policy = BrineCO2DefaultPolicy</*constantSalinity?*/true> > | ||
| 99 | class BrineCO2 | ||
| 100 | : public Base<Scalar, BrineCO2<Scalar, CO2Component, H2OType, Policy>> | ||
| 101 | , public Detail::BrineCO2Indices<Policy::useConstantSalinity()> | ||
| 102 | { | ||
| 103 | using ThisType = BrineCO2<Scalar, CO2Component, H2OType, Policy>; | ||
| 104 | // binary coefficients | ||
| 105 | using Brine_CO2 = BinaryCoeff::Brine_CO2<Scalar, CO2Component>; | ||
| 106 | |||
| 107 | // use constant salinity brine? | ||
| 108 | static constexpr bool useConstantSalinity = Policy::useConstantSalinity(); | ||
| 109 | |||
| 110 | // The possible brine types | ||
| 111 | using VariableSalinityBrine = Dumux::FluidSystems::Brine<Scalar, H2OType>; | ||
| 112 | using ConstantSalinityBrine = Dumux::Components::Brine<Scalar, H2OType>; | ||
| 113 | using BrineType = typename std::conditional_t< useConstantSalinity, | ||
| 114 | ConstantSalinityBrine, | ||
| 115 | VariableSalinityBrine >; | ||
| 116 | |||
| 117 | ///////////////////////////////////////////////////////////////////////////////// | ||
| 118 | //! The following two indices are only used internally and are not part of the | ||
| 119 | //! public interface. Depending on the chosen policy, i.e. if brine is used as | ||
| 120 | //! a pseudo component or a fluid system with NaCl as a separate component, the | ||
| 121 | //! indices that are part of the public interface are chosen by inheritance from | ||
| 122 | //! Detail::BrineCO2Indices (see documentation). | ||
| 123 | //! | ||
| 124 | //! depending on the implementation this is either brine (pseudo-component) or H2O | ||
| 125 | static constexpr int BrineOrH2OIdx = 0; | ||
| 126 | //! if the implementation considers NaCl as a real component, it gets the index 2 | ||
| 127 | static constexpr int NaClIdx = 2; | ||
| 128 | |||
| 129 | public: | ||
| 130 | using ParameterCache = NullParameterCache; | ||
| 131 | |||
| 132 | using H2O = H2OType; | ||
| 133 | using Brine = BrineType; | ||
| 134 | using CO2 = CO2Component; | ||
| 135 | |||
| 136 | static constexpr int numComponents = useConstantSalinity ? 2 : 3; | ||
| 137 | static constexpr int numPhases = 2; | ||
| 138 | |||
| 139 | static constexpr int liquidPhaseIdx = 0; //!< index of the liquid phase | ||
| 140 | static constexpr int gasPhaseIdx = 1; //!< index of the gas phase | ||
| 141 | static constexpr int phase0Idx = liquidPhaseIdx; //!< index of the first phase | ||
| 142 | static constexpr int phase1Idx = gasPhaseIdx; //!< index of the second phase | ||
| 143 | |||
| 144 | static constexpr int comp0Idx = 0; | ||
| 145 | static constexpr int comp1Idx = 1; | ||
| 146 | |||
| 147 | // CO2 is always the second component | ||
| 148 | static constexpr int CO2Idx = comp1Idx; | ||
| 149 | |||
| 150 | private: | ||
| 151 | |||
| 152 | // Adapter policy for the fluid state corresponding to the brine fluid system | ||
| 153 | struct BrineAdapterPolicy | ||
| 154 | { | ||
| 155 | using FluidSystem = VariableSalinityBrine; | ||
| 156 | |||
| 157 | ✗ | static constexpr int phaseIdx(int brinePhaseIdx) { return liquidPhaseIdx; } | |
| 158 | static constexpr int compIdx(int brineCompIdx) | ||
| 159 | { | ||
| 160 |
8/16✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 6 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 6 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 6 times.
|
32 | assert(brineCompIdx == VariableSalinityBrine::H2OIdx || brineCompIdx == VariableSalinityBrine::NaClIdx); |
| 161 |
24/32✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
✓ Branch 5 taken 6 times.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 6 times.
✓ Branch 13 taken 6 times.
✓ Branch 14 taken 6 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 2 times.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 6 times.
✓ Branch 21 taken 6 times.
✓ Branch 22 taken 6 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 2 times.
✓ Branch 25 taken 2 times.
✓ Branch 26 taken 2 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 6 times.
✓ Branch 29 taken 6 times.
✓ Branch 30 taken 6 times.
✗ Branch 31 not taken.
|
96 | switch (brineCompIdx) |
| 162 | { | ||
| 163 | case VariableSalinityBrine::H2OIdx: return BrineOrH2OIdx; | ||
| 164 | 64 | case VariableSalinityBrine::NaClIdx: return NaClIdx; | |
| 165 | default: return 0; // this will never be reached, only needed to suppress compiler warning | ||
| 166 | } | ||
| 167 | } | ||
| 168 | }; | ||
| 169 | |||
| 170 | template<class FluidState> | ||
| 171 | using BrineAdapter = FluidStateAdapter<FluidState, BrineAdapterPolicy>; | ||
| 172 | |||
| 173 | public: | ||
| 174 | |||
| 175 | /*! | ||
| 176 | * \brief Return the human readable name of a fluid phase | ||
| 177 | * \param phaseIdx The index of the fluid phase to consider | ||
| 178 | */ | ||
| 179 | 300 | static std::string phaseName(int phaseIdx) | |
| 180 | { | ||
| 181 |
2/3✓ Branch 0 taken 142 times.
✓ Branch 1 taken 142 times.
✗ Branch 2 not taken.
|
300 | switch (phaseIdx) |
| 182 | { | ||
| 183 | 150 | case liquidPhaseIdx: return IOName::liquidPhase(); | |
| 184 | 150 | case gasPhaseIdx: return IOName::gaseousPhase(); | |
| 185 | } | ||
| 186 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | |
| 187 | } | ||
| 188 | |||
| 189 | /*! | ||
| 190 | * \brief Returns whether the fluids are miscible | ||
| 191 | */ | ||
| 192 | static constexpr bool isMiscible() | ||
| 193 | { return true; } | ||
| 194 | |||
| 195 | /*! | ||
| 196 | * \brief Return whether a phase is gaseous | ||
| 197 | * | ||
| 198 | * \param phaseIdx The index of the fluid phase to consider | ||
| 199 | */ | ||
| 200 | static constexpr bool isGas(int phaseIdx) | ||
| 201 | { | ||
| 202 |
6/12✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 4 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 4 times.
|
24 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
| 203 | |||
| 204 | return phaseIdx == gasPhaseIdx; | ||
| 205 | } | ||
| 206 | |||
| 207 | /*! | ||
| 208 | * \brief Returns true if and only if a fluid phase is assumed to | ||
| 209 | * be an ideal gas. | ||
| 210 | * \param phaseIdx The index of the fluid phase to consider | ||
| 211 | */ | ||
| 212 | static constexpr bool isIdealGas(int phaseIdx) | ||
| 213 | { | ||
| 214 | assert(0 <= phaseIdx && phaseIdx < numPhases); | ||
| 215 | // let the fluids decide | ||
| 216 | if (phaseIdx == gasPhaseIdx) | ||
| 217 | return useConstantSalinity ? (ConstantSalinityBrine::gasIsIdeal() && CO2::gasIsIdeal()) | ||
| 218 | : (H2O::gasIsIdeal() && CO2::gasIsIdeal()); | ||
| 219 | return false; // not a gas | ||
| 220 | } | ||
| 221 | |||
| 222 | /*! | ||
| 223 | * \brief Returns true if and only if a fluid phase is assumed to | ||
| 224 | * be an ideal mixture. | ||
| 225 | * | ||
| 226 | * We define an ideal mixture as a fluid phase where the fugacity | ||
| 227 | * coefficients of all components times the pressure of the phase | ||
| 228 | * are independent on the fluid composition. This assumption is true | ||
| 229 | * if Henry's law and Raoult's law apply. If you are unsure what | ||
| 230 | * this function should return, it is safe to return false. The | ||
| 231 | * only damage done will be (slightly) increased computation times | ||
| 232 | * in some cases. | ||
| 233 | * | ||
| 234 | * \param phaseIdx The index of the fluid phase to consider | ||
| 235 | */ | ||
| 236 | static bool isIdealMixture(int phaseIdx) | ||
| 237 | { | ||
| 238 | assert(0 <= phaseIdx && phaseIdx < numPhases); | ||
| 239 |
16/16✓ Branch 0 taken 3 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 3 times.
✓ Branch 9 taken 3 times.
✓ Branch 10 taken 2 times.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 3 times.
✓ Branch 13 taken 3 times.
✓ Branch 14 taken 2 times.
✓ Branch 15 taken 2 times.
|
40 | if (phaseIdx == liquidPhaseIdx) |
| 240 | return false; | ||
| 241 | return true; | ||
| 242 | } | ||
| 243 | |||
| 244 | /*! | ||
| 245 | * \brief Returns true if and only if a fluid phase is assumed to | ||
| 246 | * be compressible. | ||
| 247 | * | ||
| 248 | * Compressible means that the partial derivative of the density | ||
| 249 | * to the fluid pressure is always larger than zero. | ||
| 250 | * | ||
| 251 | * \param phaseIdx The index of the fluid phase to consider | ||
| 252 | */ | ||
| 253 | static constexpr bool isCompressible(int phaseIdx) | ||
| 254 | { | ||
| 255 | assert(0 <= phaseIdx && phaseIdx < numPhases); | ||
| 256 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
|
2 | if (phaseIdx == liquidPhaseIdx) |
| 257 | return useConstantSalinity ? ConstantSalinityBrine::liquidIsCompressible() | ||
| 258 | : VariableSalinityBrine::isCompressible(VariableSalinityBrine::liquidPhaseIdx); | ||
| 259 | return true; | ||
| 260 | } | ||
| 261 | |||
| 262 | /*! | ||
| 263 | * \brief Return the human readable name of a component | ||
| 264 | * \param compIdx The index of the component to consider | ||
| 265 | */ | ||
| 266 | 152 | static std::string componentName(int compIdx) | |
| 267 | { | ||
| 268 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 132 times.
|
152 | assert(0 <= compIdx && compIdx < numComponents); |
| 269 | if (useConstantSalinity) | ||
| 270 | { | ||
| 271 |
5/10✓ Branch 0 taken 18 times.
✓ Branch 1 taken 102 times.
✓ Branch 3 taken 18 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 18 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 18 times.
✗ Branch 10 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
|
128 | static std::string name[] = { ConstantSalinityBrine::name(), CO2::name() }; |
| 272 | 128 | return name[compIdx]; | |
| 273 | } | ||
| 274 | else | ||
| 275 | { | ||
| 276 |
6/12✓ Branch 0 taken 4 times.
✓ Branch 1 taken 8 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 4 times.
✗ Branch 13 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
|
24 | static std::string name[] = { VariableSalinityBrine::componentName(VariableSalinityBrine::H2OIdx), |
| 277 | CO2::name(), | ||
| 278 | VariableSalinityBrine::NaCl::name() }; | ||
| 279 | 24 | return name[compIdx]; | |
| 280 | } | ||
| 281 | } | ||
| 282 | |||
| 283 | /*! | ||
| 284 | * \brief Return the molar mass of a component in \f$\mathrm{[kg/mol]}\f$. | ||
| 285 | * \param compIdx The index of the component to consider | ||
| 286 | */ | ||
| 287 | 5143953393 | static Scalar molarMass(int compIdx) | |
| 288 | { | ||
| 289 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5143953161 times.
|
5143953393 | assert(0 <= compIdx && compIdx < numComponents); |
| 290 | if (useConstantSalinity) | ||
| 291 | { | ||
| 292 |
4/6✓ Branch 0 taken 18 times.
✓ Branch 1 taken 5143952963 times.
✓ Branch 3 taken 18 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 18 times.
✗ Branch 7 not taken.
|
5143953033 | static const Scalar M[] = { ConstantSalinityBrine::molarMass(), CO2::molarMass() }; |
| 293 | 5143953033 | return M[compIdx]; | |
| 294 | } | ||
| 295 | else | ||
| 296 | { | ||
| 297 |
4/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 176 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
|
368 | static const Scalar M[] = { VariableSalinityBrine::molarMass(VariableSalinityBrine::H2OIdx), |
| 298 | CO2::molarMass(), | ||
| 299 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
8 | VariableSalinityBrine::molarMass(VariableSalinityBrine::NaClIdx) }; |
| 300 | 360 | return M[compIdx]; | |
| 301 | } | ||
| 302 | } | ||
| 303 | |||
| 304 | /**************************************** | ||
| 305 | * thermodynamic relations | ||
| 306 | ****************************************/ | ||
| 307 | |||
| 308 | // Initializing with salinity and default tables | ||
| 309 | static void init() | ||
| 310 | { | ||
| 311 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
10 | init(/*startTemp=*/273.15, /*endTemp=*/623.15, /*tempSteps=*/100, |
| 312 | /*startPressure=*/1e4, /*endPressure=*/40e6, /*pressureSteps=*/200); | ||
| 313 | } | ||
| 314 | |||
| 315 | // Initializing with custom tables | ||
| 316 | 32 | static void init(Scalar startTemp, Scalar endTemp, int tempSteps, | |
| 317 | Scalar startPressure, Scalar endPressure, int pressureSteps) | ||
| 318 | { | ||
| 319 | 32 | std::cout << "The Brine-CO2 fluid system was configured with the following policy:\n"; | |
| 320 | 64 | std::cout << " - use constant salinity: " << std::boolalpha << Policy::useConstantSalinity() << "\n"; | |
| 321 | 64 | std::cout << " - use CO2 gas density as gas mixture density: " << std::boolalpha << Policy::useCO2GasDensityAsGasMixtureDensity() << std::endl; | |
| 322 | |||
| 323 | if (H2O::isTabulated) | ||
| 324 | 24 | H2O::init(startTemp, endTemp, tempSteps, startPressure, endPressure, pressureSteps); | |
| 325 | 32 | } | |
| 326 | |||
| 327 | using Base<Scalar, ThisType>::density; | ||
| 328 | /*! | ||
| 329 | * \brief Given a phase's composition, temperature, pressure, and | ||
| 330 | * the partial pressures of all components, return its | ||
| 331 | * density \f$\mathrm{[kg/m^3]}\f$. | ||
| 332 | * | ||
| 333 | * \param fluidState The fluid state | ||
| 334 | * \param phaseIdx The index of the phase | ||
| 335 | */ | ||
| 336 | template <class FluidState> | ||
| 337 | 218883511 | static Scalar density(const FluidState& fluidState, int phaseIdx) | |
| 338 | { | ||
| 339 |
2/2✓ Branch 0 taken 145035962 times.
✓ Branch 1 taken 73847501 times.
|
218883511 | Scalar T = fluidState.temperature(phaseIdx); |
| 340 |
2/2✓ Branch 0 taken 145035978 times.
✓ Branch 1 taken 73847509 times.
|
218883511 | if (phaseIdx == liquidPhaseIdx) |
| 341 | 145035994 | return liquidDensityMixture_(fluidState); | |
| 342 | |||
| 343 |
1/2✓ Branch 0 taken 73847509 times.
✗ Branch 1 not taken.
|
73847517 | else if (phaseIdx == gasPhaseIdx) |
| 344 | { | ||
| 345 | if (Policy::useCO2GasDensityAsGasMixtureDensity()) | ||
| 346 | // use the CO2 gas density only and neglect compositional effects | ||
| 347 | 213565387 | return CO2::gasDensity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); | |
| 348 | else | ||
| 349 | { | ||
| 350 | // assume ideal mixture: steam and CO2 don't "see" each other | ||
| 351 | 2659052 | Scalar rho_gH2O = H2O::gasDensity(T, fluidState.partialPressure(gasPhaseIdx, BrineOrH2OIdx)); | |
| 352 | 5318092 | Scalar rho_gCO2 = CO2::gasDensity(T, fluidState.partialPressure(gasPhaseIdx, CO2Idx)); | |
| 353 | 2659052 | return (rho_gH2O + rho_gCO2); | |
| 354 | } | ||
| 355 | } | ||
| 356 | |||
| 357 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index."); | |
| 358 | } | ||
| 359 | |||
| 360 | using Base<Scalar, ThisType>::molarDensity; | ||
| 361 | //! \copydoc Base<Scalar,ThisType>::molarDensity(const FluidState&,int) | ||
| 362 | template <class FluidState> | ||
| 363 | 142376954 | static Scalar molarDensity(const FluidState& fluidState, int phaseIdx) | |
| 364 | { | ||
| 365 |
2/2✓ Branch 0 taken 71188461 times.
✓ Branch 1 taken 71188461 times.
|
142376954 | Scalar T = fluidState.temperature(phaseIdx); |
| 366 |
2/2✓ Branch 0 taken 71188469 times.
✓ Branch 1 taken 71188469 times.
|
142376954 | if (phaseIdx == liquidPhaseIdx) |
| 367 | 71188477 | return density(fluidState, phaseIdx)/fluidState.averageMolarMass(phaseIdx); | |
| 368 |
1/2✓ Branch 0 taken 71188469 times.
✗ Branch 1 not taken.
|
71188477 | else if (phaseIdx == gasPhaseIdx) |
| 369 | { | ||
| 370 | if (Policy::useCO2GasDensityAsGasMixtureDensity()) | ||
| 371 | 213565387 | return CO2::gasMolarDensity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); | |
| 372 | else | ||
| 373 | { | ||
| 374 | // assume ideal mixture: steam and CO2 don't "see" each other | ||
| 375 | 12 | Scalar rhoMolar_gH2O = H2O::gasMolarDensity(T, fluidState.partialPressure(gasPhaseIdx, BrineOrH2OIdx)); | |
| 376 | 12 | Scalar rhoMolar_gCO2 = CO2::gasMolarDensity(T, fluidState.partialPressure(gasPhaseIdx, CO2Idx)); | |
| 377 | 12 | return rhoMolar_gH2O + rhoMolar_gCO2; | |
| 378 | } | ||
| 379 | } | ||
| 380 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index."); | |
| 381 | } | ||
| 382 | |||
| 383 | using Base<Scalar, ThisType>::viscosity; | ||
| 384 | /*! | ||
| 385 | * \brief Calculate the dynamic viscosity of a fluid phase \f$\mathrm{[Pa*s]}\f$ | ||
| 386 | * | ||
| 387 | * \param fluidState An arbitrary fluid state | ||
| 388 | * \param phaseIdx The index of the fluid phase to consider | ||
| 389 | * | ||
| 390 | * \note For the viscosity of the phases the contribution of the minor | ||
| 391 | * component is neglected. This contribution is probably not big, but somebody | ||
| 392 | * would have to find out its influence. | ||
| 393 | */ | ||
| 394 | template <class FluidState> | ||
| 395 | 147695034 | static Scalar viscosity(const FluidState& fluidState, int phaseIdx) | |
| 396 | { | ||
| 397 |
2/2✓ Branch 0 taken 73847501 times.
✓ Branch 1 taken 73847501 times.
|
147695034 | Scalar T = fluidState.temperature(phaseIdx); |
| 398 |
2/2✓ Branch 0 taken 73847501 times.
✓ Branch 1 taken 73847501 times.
|
147695034 | Scalar p = fluidState.pressure(phaseIdx); |
| 399 | |||
| 400 |
2/2✓ Branch 0 taken 73847509 times.
✓ Branch 1 taken 73847509 times.
|
147695034 | if (phaseIdx == liquidPhaseIdx) |
| 401 | 73847525 | return useConstantSalinity ? ConstantSalinityBrine::liquidViscosity(T, p) | |
| 402 | : VariableSalinityBrine::viscosity( BrineAdapter<FluidState>(fluidState), | ||
| 403 | 73847509 | VariableSalinityBrine::liquidPhaseIdx ); | |
| 404 |
1/2✓ Branch 0 taken 73847509 times.
✗ Branch 1 not taken.
|
73847517 | else if (phaseIdx == gasPhaseIdx) |
| 405 | 73847517 | return CO2::gasViscosity(T, p); | |
| 406 | |||
| 407 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index."); | |
| 408 | } | ||
| 409 | |||
| 410 | using Base<Scalar, ThisType>::fugacityCoefficient; | ||
| 411 | //! \copydoc Base<Scalar,ThisType>::fugacityCoefficient(const FluidState&,int,int) | ||
| 412 | template <class FluidState> | ||
| 413 | 80 | static Scalar fugacityCoefficient(const FluidState& fluidState, | |
| 414 | int phaseIdx, | ||
| 415 | int compIdx) | ||
| 416 | { | ||
| 417 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
|
80 | assert(0 <= compIdx && compIdx < numComponents); |
| 418 | |||
| 419 |
2/2✓ Branch 0 taken 20 times.
✓ Branch 1 taken 20 times.
|
80 | if (phaseIdx == gasPhaseIdx) |
| 420 | // use the fugacity coefficients of an ideal gas. the | ||
| 421 | // actual value of the fugacity is not relevant, as long | ||
| 422 | // as the relative fluid compositions are observed, | ||
| 423 | return 1.0; | ||
| 424 | |||
| 425 |
1/2✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
|
40 | else if (phaseIdx == liquidPhaseIdx) |
| 426 | { | ||
| 427 | 40 | Scalar T = fluidState.temperature(phaseIdx); | |
| 428 | 40 | Scalar pl = fluidState.pressure(liquidPhaseIdx); | |
| 429 | 40 | Scalar pg = fluidState.pressure(gasPhaseIdx); | |
| 430 | |||
| 431 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 20 times.
|
40 | assert(T > 0); |
| 432 |
2/4✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 20 times.
|
40 | assert(pl > 0 && pg > 0); |
| 433 | |||
| 434 | // calculate the equilibrium composition for given T & p | ||
| 435 | Scalar xlH2O, xgH2O; | ||
| 436 | Scalar xlCO2, xgCO2; | ||
| 437 | 40 | const Scalar salinity = useConstantSalinity ? ConstantSalinityBrine::salinity() | |
| 438 | : fluidState.massFraction(liquidPhaseIdx, NaClIdx); | ||
| 439 | 40 | Brine_CO2::calculateMoleFractions(T, pl, salinity, /*knownGasPhaseIdx=*/-1, xlCO2, xgH2O); | |
| 440 | |||
| 441 | // normalize the phase compositions | ||
| 442 | using std::min; | ||
| 443 | using std::max; | ||
| 444 |
2/4✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
|
80 | xlCO2 = max(0.0, min(1.0, xlCO2)); |
| 445 |
2/4✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
|
80 | xgH2O = max(0.0, min(1.0, xgH2O)); |
| 446 | 40 | xlH2O = 1.0 - xlCO2; | |
| 447 | 40 | xgCO2 = 1.0 - xgH2O; | |
| 448 | |||
| 449 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 12 times.
|
40 | if (compIdx == BrineOrH2OIdx) |
| 450 | 16 | return (xgH2O/xlH2O)*(pg/pl); | |
| 451 | |||
| 452 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
|
16 | else if (compIdx == CO2Idx) |
| 453 | 16 | return (xgCO2/xlCO2)*(pg/pl); | |
| 454 | |||
| 455 | // NaCl is assumed to stay in the liquid! | ||
| 456 | else if (!useConstantSalinity && compIdx == NaClIdx) | ||
| 457 | return 0.0; | ||
| 458 | |||
| 459 | DUNE_THROW(Dune::InvalidStateException, "Invalid component index."); | ||
| 460 | } | ||
| 461 | |||
| 462 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index."); | |
| 463 | } | ||
| 464 | |||
| 465 | /*! | ||
| 466 | * \brief Returns the equilibrium mole fraction of the dissolved component in a phase. | ||
| 467 | * \param fluidState An arbitrary fluid state | ||
| 468 | * \param paramCache Parameter cache | ||
| 469 | * \param phaseIdx The index of the fluid phase to consider | ||
| 470 | */ | ||
| 471 | template <class FluidState> | ||
| 472 | 80456217 | static Scalar equilibriumMoleFraction(const FluidState& fluidState, | |
| 473 | const ParameterCache& paramCache, | ||
| 474 | int phaseIdx) | ||
| 475 | { | ||
| 476 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 80456217 times.
|
80456217 | Scalar T = fluidState.temperature(phaseIdx); |
| 477 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 80456217 times.
|
80456217 | Scalar p = fluidState.pressure(phaseIdx); |
| 478 | |||
| 479 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 80456217 times.
|
80456217 | assert(T > 0); |
| 480 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 80456217 times.
|
80456217 | assert(p > 0); |
| 481 | |||
| 482 | Scalar xgH2O; | ||
| 483 | Scalar xlCO2; | ||
| 484 | |||
| 485 | // calculate the equilibrium composition for given T & p | ||
| 486 | 80456217 | const Scalar salinity = useConstantSalinity ? ConstantSalinityBrine::salinity() | |
| 487 | : fluidState.massFraction(liquidPhaseIdx, NaClIdx); | ||
| 488 | 80456217 | Brine_CO2::calculateMoleFractions(T, p, salinity, /*knowgasPhaseIdx=*/-1, xlCO2, xgH2O); | |
| 489 | |||
| 490 |
2/2✓ Branch 0 taken 71188461 times.
✓ Branch 1 taken 9267756 times.
|
80456217 | if (phaseIdx == gasPhaseIdx) |
| 491 | 71188461 | return xgH2O; | |
| 492 |
1/2✓ Branch 0 taken 9267756 times.
✗ Branch 1 not taken.
|
9267756 | else if (phaseIdx == liquidPhaseIdx) |
| 493 | 9267756 | return xlCO2; | |
| 494 | |||
| 495 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index."); | |
| 496 | } | ||
| 497 | |||
| 498 | |||
| 499 | using Base<Scalar, ThisType>::diffusionCoefficient; | ||
| 500 | /*! | ||
| 501 | * \brief Calculate the molecular diffusion coefficient for a | ||
| 502 | * component in a fluid phase \f$\mathrm{[mol^2 * s / (kg*m^3)]}\f$ | ||
| 503 | * | ||
| 504 | * Molecular diffusion of a component \f$\mathrm{\kappa}\f$ is caused by a | ||
| 505 | * gradient of the chemical potential and follows the law | ||
| 506 | * | ||
| 507 | * \f[ J = - D \nabla mu_\kappa \f] | ||
| 508 | * | ||
| 509 | * where \f$\mathrm{\mu_\kappa}\f$ is the component's chemical potential, | ||
| 510 | * \f$D\f$ is the diffusion coefficient and \f$\mathrm{J}\f$ is the | ||
| 511 | * diffusive flux. \f$\mathrm{mu_\kappa}\f$ is connected to the component's | ||
| 512 | * fugacity \f$\mathrm{f_\kappa}\f$ by the relation | ||
| 513 | * | ||
| 514 | * \f[ \mu_\kappa = R T_\alpha \mathrm{ln} \frac{f_\kappa}{p_\alpha} \f] | ||
| 515 | * | ||
| 516 | * where \f$\mathrm{p_\alpha}\f$ and \f$\mathrm{T_\alpha}\f$ are the fluid phase' | ||
| 517 | * pressure and temperature. | ||
| 518 | * | ||
| 519 | * Maybe see http://www.ddbst.de/en/EED/PCP/DIF_C1050.php | ||
| 520 | * | ||
| 521 | * \param fluidState An arbitrary fluid state | ||
| 522 | * \param phaseIdx The index of the fluid phase to consider | ||
| 523 | * \param compIdx The index of the component to consider | ||
| 524 | */ | ||
| 525 | template <class FluidState> | ||
| 526 | 80 | static Scalar diffusionCoefficient(const FluidState& fluidState, int phaseIdx, int compIdx) | |
| 527 |
7/16✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 40 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 40 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 40 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 40 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 40 times.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 40 times.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
|
880 | { DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients"); } |
| 528 | |||
| 529 | using Base<Scalar, ThisType>::binaryDiffusionCoefficient; | ||
| 530 | // \copydoc Base<Scalar,ThisType>::binaryDiffusionCoefficient(const FluidState&,int,int,int) | ||
| 531 | template <class FluidState> | ||
| 532 | 142377130 | static Scalar binaryDiffusionCoefficient(const FluidState& fluidState, | |
| 533 | int phaseIdx, | ||
| 534 | int compIIdx, | ||
| 535 | int compJIdx) | ||
| 536 | { | ||
| 537 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 142377026 times.
|
142377130 | assert(0 <= compIIdx && compIIdx < numComponents); |
| 538 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 142377026 times.
|
142377130 | assert(0 <= compJIdx && compJIdx < numComponents); |
| 539 | |||
| 540 |
2/2✓ Branch 0 taken 71188493 times.
✓ Branch 1 taken 71188533 times.
|
142377130 | if (compIIdx > compJIdx) |
| 541 | { | ||
| 542 | using std::swap; | ||
| 543 | 71188525 | swap(compIIdx, compJIdx); | |
| 544 | } | ||
| 545 | |||
| 546 |
2/2✓ Branch 0 taken 71188461 times.
✓ Branch 1 taken 71188461 times.
|
142377130 | Scalar T = fluidState.temperature(phaseIdx); |
| 547 |
2/2✓ Branch 0 taken 71188461 times.
✓ Branch 1 taken 71188461 times.
|
142377130 | Scalar p = fluidState.pressure(phaseIdx); |
| 548 |
2/2✓ Branch 0 taken 71188513 times.
✓ Branch 1 taken 71188513 times.
|
142377130 | if (phaseIdx == liquidPhaseIdx) |
| 549 | { | ||
| 550 |
4/4✓ Branch 0 taken 71188493 times.
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 71188477 times.
✓ Branch 3 taken 16 times.
|
71188565 | if (compIIdx == BrineOrH2OIdx && compJIdx == CO2Idx) |
| 551 | 71188493 | return Brine_CO2::liquidDiffCoeff(T, p); | |
| 552 |
4/4✓ Branch 0 taken 12 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 4 times.
|
56 | if (!useConstantSalinity && compIIdx == BrineOrH2OIdx && compJIdx == NaClIdx) |
| 553 | 32 | return VariableSalinityBrine::binaryDiffusionCoefficient( BrineAdapter<FluidState>(fluidState), | |
| 554 | VariableSalinityBrine::liquidPhaseIdx, | ||
| 555 | VariableSalinityBrine::H2OIdx, | ||
| 556 | VariableSalinityBrine::NaClIdx ); | ||
| 557 | |||
| 558 |
9/20✓ Branch 2 taken 28 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 28 times.
✗ Branch 12 not taken.
✓ Branch 16 taken 28 times.
✗ Branch 17 not taken.
✓ Branch 20 taken 28 times.
✗ Branch 21 not taken.
✓ Branch 24 taken 28 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 28 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 28 times.
✗ Branch 31 not taken.
✓ Branch 32 taken 28 times.
✗ Branch 33 not taken.
✗ Branch 35 not taken.
✓ Branch 36 taken 28 times.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
|
784 | DUNE_THROW(Dune::NotImplemented, "Binary diffusion coefficient of components " << |
| 559 | compIIdx << " and " << compJIdx << " in phase " << phaseIdx); | ||
| 560 | } | ||
| 561 |
1/2✓ Branch 0 taken 71188513 times.
✗ Branch 1 not taken.
|
71188565 | else if (phaseIdx == gasPhaseIdx) |
| 562 | { | ||
| 563 |
4/4✓ Branch 0 taken 71188493 times.
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 71188477 times.
✓ Branch 3 taken 16 times.
|
71188565 | if (compIIdx == BrineOrH2OIdx && compJIdx == CO2Idx) |
| 564 | 71188493 | return Brine_CO2::gasDiffCoeff(T, p); | |
| 565 | |||
| 566 | // NaCl is expected to never be present in the gas phase. we need to | ||
| 567 | // return a diffusion coefficient that does not case numerical problems. | ||
| 568 | // We choose a very small value here. | ||
| 569 |
4/4✓ Branch 0 taken 12 times.
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 8 times.
|
56 | else if (!useConstantSalinity && compIIdx == CO2Idx && compJIdx == NaClIdx) |
| 570 | return 1e-12; | ||
| 571 | |||
| 572 |
9/20✓ Branch 2 taken 28 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 28 times.
✗ Branch 12 not taken.
✓ Branch 16 taken 28 times.
✗ Branch 17 not taken.
✓ Branch 20 taken 28 times.
✗ Branch 21 not taken.
✓ Branch 24 taken 28 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 28 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 28 times.
✗ Branch 31 not taken.
✓ Branch 32 taken 28 times.
✗ Branch 33 not taken.
✗ Branch 35 not taken.
✓ Branch 36 taken 28 times.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
|
784 | DUNE_THROW(Dune::NotImplemented, "Binary diffusion coefficient of components " << |
| 573 | compIIdx << " and " << compJIdx << " in phase " << phaseIdx); | ||
| 574 | } | ||
| 575 | |||
| 576 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index."); | |
| 577 | } | ||
| 578 | |||
| 579 | using Base<Scalar, ThisType>::enthalpy; | ||
| 580 | // \copydoc Base<Scalar,ThisType>::enthalpy(const FluidState&,int) | ||
| 581 | template <class FluidState> | ||
| 582 | 92687060 | static Scalar enthalpy(const FluidState& fluidState, int phaseIdx) | |
| 583 | { | ||
| 584 |
2/2✓ Branch 0 taken 46343514 times.
✓ Branch 1 taken 46343514 times.
|
92687060 | Scalar T = fluidState.temperature(phaseIdx); |
| 585 |
2/2✓ Branch 0 taken 46343514 times.
✓ Branch 1 taken 46343514 times.
|
92687060 | Scalar p = fluidState.pressure(phaseIdx); |
| 586 | |||
| 587 |
2/2✓ Branch 0 taken 46343522 times.
✓ Branch 1 taken 46343522 times.
|
92687060 | if (phaseIdx == liquidPhaseIdx) |
| 588 | { | ||
| 589 | // Convert J/kg to kJ/kg | ||
| 590 | 46343530 | const Scalar h_ls1 = useConstantSalinity ? ConstantSalinityBrine::liquidEnthalpy(T, p)/1e3 | |
| 591 | 16 | : VariableSalinityBrine::enthalpy( BrineAdapter<FluidState>(fluidState), | |
| 592 | VariableSalinityBrine::liquidPhaseIdx )/1e3; | ||
| 593 | |||
| 594 | // mass fraction of CO2 in Brine | ||
| 595 | 46343530 | const Scalar X_CO2_w = fluidState.massFraction(liquidPhaseIdx, CO2Idx); | |
| 596 | |||
| 597 | // heat of dissolution for CO2 according to Fig. 6 in Duan and Sun 2003. (kJ/kg) | ||
| 598 | // In the relevant temperature ranges CO2 dissolution is exothermal | ||
| 599 | 46343530 | const Scalar delta_hCO2 = (-57.4375 + T * 0.1325) * 1000/44; | |
| 600 | |||
| 601 | // enthalpy contribution of water and CO2 (kJ/kg) | ||
| 602 | 46343530 | const Scalar hw = H2O::liquidEnthalpy(T, p)/1e3; | |
| 603 | 46343530 | const Scalar hg = CO2::liquidEnthalpy(T, p)/1e3 + delta_hCO2; | |
| 604 | |||
| 605 | // Enthalpy of brine with dissolved CO2 (kJ/kg) | ||
| 606 | 46343530 | return (h_ls1 - X_CO2_w*hw + hg*X_CO2_w)*1e3; | |
| 607 | } | ||
| 608 |
1/2✓ Branch 0 taken 46343522 times.
✗ Branch 1 not taken.
|
46343530 | else if (phaseIdx == gasPhaseIdx) |
| 609 | { | ||
| 610 | // we assume NaCl to not enter the gas phase, only consider H2O and CO2 | ||
| 611 | 46343530 | return H2O::gasEnthalpy(T, p)*fluidState.massFraction(gasPhaseIdx, BrineOrH2OIdx) | |
| 612 | 46343530 | + CO2::gasEnthalpy(T, p) *fluidState.massFraction(gasPhaseIdx, CO2Idx); | |
| 613 | } | ||
| 614 | |||
| 615 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index."); | |
| 616 | } | ||
| 617 | |||
| 618 | /*! | ||
| 619 | * \brief Returns the specific enthalpy \f$\mathrm{[J/kg]}\f$ of a component in a specific phase | ||
| 620 | * \param fluidState The fluid state | ||
| 621 | * \param phaseIdx The index of the phase | ||
| 622 | * \param componentIdx The index of the component | ||
| 623 | */ | ||
| 624 | template <class FluidState> | ||
| 625 | static Scalar componentEnthalpy(const FluidState& fluidState, int phaseIdx, int componentIdx) | ||
| 626 | { | ||
| 627 | const Scalar T = fluidState.temperature(phaseIdx); | ||
| 628 | const Scalar p = fluidState.pressure(phaseIdx); | ||
| 629 | |||
| 630 | if (phaseIdx == liquidPhaseIdx) | ||
| 631 | { | ||
| 632 | if (componentIdx == BrineOrH2OIdx) | ||
| 633 | return H2O::liquidEnthalpy(T, p); | ||
| 634 | else if (componentIdx == CO2Idx) | ||
| 635 | return CO2::liquidEnthalpy(T, p); | ||
| 636 | else if (componentIdx == NaClIdx) | ||
| 637 | DUNE_THROW(Dune::NotImplemented, "The component enthalpy for NaCl is not implemented."); | ||
| 638 | DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << componentIdx); | ||
| 639 | } | ||
| 640 | else if (phaseIdx == gasPhaseIdx) | ||
| 641 | { | ||
| 642 | if (componentIdx == BrineOrH2OIdx) | ||
| 643 | return H2O::gasEnthalpy(T, p); | ||
| 644 | else if (componentIdx == CO2Idx) | ||
| 645 | return CO2::gasEnthalpy(T, p); | ||
| 646 | else if (componentIdx == NaClIdx) | ||
| 647 | DUNE_THROW(Dune::InvalidStateException, "Implementation assumes NaCl not to be present in gas phase"); | ||
| 648 | DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << componentIdx); | ||
| 649 | } | ||
| 650 | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | ||
| 651 | } | ||
| 652 | |||
| 653 | using Base<Scalar, ThisType>::thermalConductivity; | ||
| 654 | /*! | ||
| 655 | * \brief Thermal conductivity of a fluid phase \f$\mathrm{[W/(m K)]}\f$. | ||
| 656 | * \param fluidState An arbitrary fluid state | ||
| 657 | * \param phaseIdx The index of the fluid phase to consider | ||
| 658 | * | ||
| 659 | * \note For the thermal conductivity of the phases the contribution of the minor | ||
| 660 | * component is neglected. This contribution is probably not big, but somebody | ||
| 661 | * would have to find out its influence. | ||
| 662 | */ | ||
| 663 | template <class FluidState> | ||
| 664 | 92687060 | static Scalar thermalConductivity(const FluidState& fluidState, int phaseIdx) | |
| 665 | { | ||
| 666 |
2/2✓ Branch 0 taken 46343522 times.
✓ Branch 1 taken 46343522 times.
|
92687060 | if (phaseIdx == liquidPhaseIdx) |
| 667 | 139030566 | return useConstantSalinity ? ConstantSalinityBrine::liquidThermalConductivity( fluidState.temperature(phaseIdx), | |
| 668 | fluidState.pressure(phaseIdx) ) | ||
| 669 | : VariableSalinityBrine::thermalConductivity( BrineAdapter<FluidState>(fluidState), | ||
| 670 | 46343522 | VariableSalinityBrine::liquidPhaseIdx ); | |
| 671 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 46343514 times.
|
46343530 | else if (phaseIdx == gasPhaseIdx) |
| 672 | 16 | return CO2::gasThermalConductivity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); | |
| 673 | |||
| 674 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index."); | |
| 675 | } | ||
| 676 | |||
| 677 | using Base<Scalar, ThisType>::heatCapacity; | ||
| 678 | /*! | ||
| 679 | * \copybrief Base<Scalar,ThisType>::heatCapacity(const FluidState&,int) | ||
| 680 | * | ||
| 681 | * \note We employ the heat capacity of the pure phases. | ||
| 682 | * | ||
| 683 | * \todo TODO Implement heat capacity for gaseous CO2 | ||
| 684 | * | ||
| 685 | * \param fluidState An arbitrary fluid state | ||
| 686 | * \param phaseIdx The index of the fluid phase to consider | ||
| 687 | */ | ||
| 688 | template <class FluidState> | ||
| 689 | 32 | static Scalar heatCapacity(const FluidState &fluidState, | |
| 690 | int phaseIdx) | ||
| 691 | { | ||
| 692 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
|
32 | if(phaseIdx == liquidPhaseIdx) |
| 693 | 24 | return useConstantSalinity ? ConstantSalinityBrine::liquidHeatCapacity( fluidState.temperature(phaseIdx), | |
| 694 | fluidState.pressure(phaseIdx) ) | ||
| 695 | : VariableSalinityBrine::heatCapacity( BrineAdapter<FluidState>(fluidState), | ||
| 696 | 8 | VariableSalinityBrine::liquidPhaseIdx ); | |
| 697 |
1/2✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
|
16 | else if (phaseIdx == gasPhaseIdx) |
| 698 | 16 | return CO2::liquidHeatCapacity(fluidState.temperature(phaseIdx), | |
| 699 | 16 | fluidState.pressure(phaseIdx)); | |
| 700 | |||
| 701 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index."); | |
| 702 | } | ||
| 703 | |||
| 704 | private: | ||
| 705 | |||
| 706 | /*! | ||
| 707 | * \brief Liquid-phase density calculation of a mixture of brine and CO2 accounting for compositional effects. | ||
| 708 | * | ||
| 709 | * \param fluidState An arbitrary fluid state | ||
| 710 | * \return liquidDensity the liquid-phase density | ||
| 711 | */ | ||
| 712 | template<class FluidState> | ||
| 713 | 145035994 | static Scalar liquidDensityMixture_(const FluidState& fluidState) | |
| 714 | { | ||
| 715 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 145035962 times.
|
145035994 | const auto T = fluidState.temperature(liquidPhaseIdx); |
| 716 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 145035962 times.
|
145035994 | const auto p = fluidState.pressure(liquidPhaseIdx); |
| 717 | |||
| 718 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 145035978 times.
|
145035994 | if (T < 273.15) |
| 719 | ✗ | DUNE_THROW(NumericalProblem, "Liquid density for Brine and CO2 is only " | |
| 720 | "defined above 273.15K (T = " << T << ")"); | ||
| 721 | |||
| 722 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 145035978 times.
|
145035994 | if (p >= 2.5e8) |
| 723 | ✗ | DUNE_THROW(NumericalProblem, "Liquid density for Brine and CO2 is only " | |
| 724 | "defined below 250MPa (p = " << p << ")"); | ||
| 725 | |||
| 726 | // density of pure water | ||
| 727 | 145035994 | Scalar rho_pure = H2O::liquidDensity(T, p); | |
| 728 | |||
| 729 | // density of water with dissolved CO2 (neglect NaCl) | ||
| 730 | Scalar rho_lCO2; | ||
| 731 | if (useConstantSalinity) | ||
| 732 | { | ||
| 733 | // use normalized composition for to calculate the density | ||
| 734 | // (the relations don't seem to take non-normalized compositions too well...) | ||
| 735 | // TODO Do we really need this normalization??? | ||
| 736 | using std::min; | ||
| 737 | using std::max; | ||
| 738 |
6/6✓ Branch 0 taken 142376922 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 142376922 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 51624166 times.
✓ Branch 5 taken 90752756 times.
|
429789838 | Scalar xlBrine = min(1.0, max(0.0, fluidState.moleFraction(liquidPhaseIdx, BrineOrH2OIdx))); |
| 739 |
5/6✓ Branch 0 taken 117565108 times.
✓ Branch 1 taken 24811822 times.
✓ Branch 2 taken 117565108 times.
✓ Branch 3 taken 24811822 times.
✓ Branch 4 taken 142376922 times.
✗ Branch 5 not taken.
|
404978024 | Scalar xlCO2 = min(1.0, max(0.0, fluidState.moleFraction(liquidPhaseIdx, CO2Idx))); |
| 740 | 145035978 | Scalar sumx = xlBrine + xlCO2; | |
| 741 | 145035978 | xlBrine /= sumx; | |
| 742 | 145035978 | xlCO2 /= sumx; | |
| 743 | |||
| 744 | 145035982 | rho_lCO2 = liquidDensityWaterCO2_(T, p, xlBrine, xlCO2); | |
| 745 | } | ||
| 746 | else | ||
| 747 | { | ||
| 748 | // rescale mole fractions | ||
| 749 | 16 | auto xlH2O = fluidState.moleFraction(liquidPhaseIdx, BrineOrH2OIdx); | |
| 750 | 16 | auto xlCO2 = fluidState.moleFraction(liquidPhaseIdx, NaClIdx); | |
| 751 | 16 | const auto sumMoleFrac = xlH2O + xlCO2; | |
| 752 | 16 | xlH2O = xlH2O/sumMoleFrac; | |
| 753 | 16 | xlCO2 = xlCO2/sumMoleFrac; | |
| 754 | 20 | rho_lCO2 = liquidDensityWaterCO2_(T, p, xlH2O, xlCO2); | |
| 755 | } | ||
| 756 | |||
| 757 | // density of brine (water with nacl) | ||
| 758 | 145036010 | Scalar rho_brine = useConstantSalinity ? ConstantSalinityBrine::liquidDensity(T, p) | |
| 759 | : VariableSalinityBrine::density( BrineAdapter<FluidState>(fluidState), | ||
| 760 | VariableSalinityBrine::liquidPhaseIdx ); | ||
| 761 | |||
| 762 | // contribution of co2 to the density | ||
| 763 | 145035994 | Scalar contribCO2 = rho_lCO2 - rho_pure; | |
| 764 | |||
| 765 | // Total brine density with dissolved CO2 | ||
| 766 | // rho_{b, CO2} = rho_pure + contribution(salt) + contribution(CO2) | ||
| 767 | 145035994 | return rho_brine + contribCO2; | |
| 768 | } | ||
| 769 | |||
| 770 | |||
| 771 | /*! | ||
| 772 | * \brief Liquid-phase density for a mixture of CO2 in pure water. | ||
| 773 | * \note this is used by liquidDensityMixture_ | ||
| 774 | * | ||
| 775 | * \param temperature The temperature | ||
| 776 | * \param pl the liquid-phase pressure | ||
| 777 | * \param xlH2O the liquid-phase H2O mole fraction | ||
| 778 | * \param xlCO2 the liquid-phase CO2 mole fraction | ||
| 779 | * \return the density of a mixture of CO2 in pure water | ||
| 780 | */ | ||
| 781 | ✗ | static Scalar liquidDensityWaterCO2_(Scalar temperature, | |
| 782 | Scalar pl, | ||
| 783 | Scalar xlH2O, | ||
| 784 | Scalar xlCO2) | ||
| 785 | { | ||
| 786 | 4 | const Scalar M_CO2 = CO2::molarMass(); | |
| 787 | 4 | const Scalar M_H2O = H2O::molarMass(); | |
| 788 | |||
| 789 | 4 | const Scalar tempC = temperature - 273.15; /* tempC : temperature in °C */ | |
| 790 | 4 | const Scalar rho_pure = H2O::liquidDensity(temperature, pl); | |
| 791 | |||
| 792 | // xlH2O is available, but in case of a pure gas phase | ||
| 793 | // the value of M_T for the virtual liquid phase can become very large | ||
| 794 | 4 | xlH2O = 1.0 - xlCO2; | |
| 795 | 4 | const Scalar M_T = M_H2O * xlH2O + M_CO2 * xlCO2; | |
| 796 | 4 | const Scalar V_phi = (37.51 + | |
| 797 | 4 | tempC*(-9.585e-2 + | |
| 798 | 4 | tempC*(8.74e-4 - | |
| 799 | 4 | tempC*5.044e-7))) / 1.0e6; | |
| 800 | 4 | return 1/(xlCO2 * V_phi/M_T + M_H2O * xlH2O / (rho_pure * M_T)); | |
| 801 | } | ||
| 802 | }; | ||
| 803 | |||
| 804 | } // end namespace Dumux::FluidSystems | ||
| 805 | |||
| 806 | #endif | ||
| 807 |