| 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 | #ifndef DUMUX_BIOFLUID_SIMPLE_CHEM_SYSTEM_HH | ||
| 9 | #define DUMUX_BIOFLUID_SIMPLE_CHEM_SYSTEM_HH | ||
| 10 | |||
| 11 | // ## The fluid system (`biominsimplechemistry.hh`) | ||
| 12 | // | ||
| 13 | // This file contains the __fluidsystem class__ which defines all functions needed to describe the fluids and their properties, | ||
| 14 | // which are needed to model biomineralization. | ||
| 15 | // | ||
| 16 | // [[content]] | ||
| 17 | // | ||
| 18 | // ### Include files | ||
| 19 | // [[details]] | ||
| 20 | // [[codeblock]] | ||
| 21 | #include <dumux/material/idealgas.hh> | ||
| 22 | |||
| 23 | // we include the base fluid system | ||
| 24 | #include <dumux/material/fluidsystems/base.hh> | ||
| 25 | // we include the brine adapter adapting component indices to use the brine fluidsystem | ||
| 26 | #include "icpcomplexsalinitybrine.hh" | ||
| 27 | |||
| 28 | // we include all necessary fluid components | ||
| 29 | #include <dumux/material/fluidstates/adapter.hh> | ||
| 30 | #include <dumux/material/components/simpleco2.hh> | ||
| 31 | #include <dumux/material/components/h2o.hh> | ||
| 32 | #include <dumux/material/components/tabulatedcomponent.hh> | ||
| 33 | #include <dumux/material/components/sodiumion.hh> | ||
| 34 | #include <dumux/material/components/chlorideion.hh> | ||
| 35 | #include <dumux/material/components/calciumion.hh> | ||
| 36 | #include <dumux/material/components/urea.hh> | ||
| 37 | #include <dumux/material/components/o2.hh> | ||
| 38 | #include <dumux/material/components/glucose.hh> | ||
| 39 | #include <examples/biomineralization/material/components/suspendedbiomass.hh> | ||
| 40 | |||
| 41 | // we include brine-co2 and h2o-co2 binary coefficients | ||
| 42 | #include <dumux/material/binarycoefficients/brine_co2.hh> | ||
| 43 | #include <dumux/material/binarycoefficients/h2o_o2.hh> | ||
| 44 | // [[/codeblock]] | ||
| 45 | |||
| 46 | #include <dumux/material/fluidsystems/nullparametercache.hh> | ||
| 47 | #include <dumux/common/exceptions.hh> | ||
| 48 | |||
| 49 | #include <assert.h> | ||
| 50 | // [[/details]] | ||
| 51 | // | ||
| 52 | // ### The fluidsystem class | ||
| 53 | // In the BioMinSimpleChemistryFluid fluid system, we define all functions needed to describe the fluids and their properties accounted for in our simulation. | ||
| 54 | // The simplified biogeochemistry biomineralization fluid system requires the CO2 component and the H2OType as template parameters. | ||
| 55 | // We enter the namespace Dumux. All Dumux functions and classes are in a namespace Dumux, to make sure they don`t clash with symbols from other libraries you may want to use in conjunction with Dumux. | ||
| 56 | // [[codeblock]] | ||
| 57 | namespace Dumux::FluidSystems { | ||
| 58 | |||
| 59 | template <class Scalar, | ||
| 60 | class CO2Impl = Components::SimpleCO2<Scalar>, | ||
| 61 | class H2OType = Components::TabulatedComponent<Components::H2O<Scalar>> > | ||
| 62 | class BioMinSimpleChemistryFluid | ||
| 63 | : public Base<Scalar, BioMinSimpleChemistryFluid<Scalar, CO2Impl, H2OType> > | ||
| 64 | { | ||
| 65 | using ThisType = BioMinSimpleChemistryFluid<Scalar, CO2Impl, H2OType>; | ||
| 66 | using Base = Dumux::FluidSystems::Base<Scalar, ThisType>; | ||
| 67 | using IdealGas = Dumux::IdealGas<Scalar>; | ||
| 68 | // [[/codeblock]] | ||
| 69 | // | ||
| 70 | // #### Component and phase definitions | ||
| 71 | // With the following function we define what phases and components will be used by the fluid system and define the indices used to distinguish phases and components in the course of the simulation. | ||
| 72 | // [[codeblock]] | ||
| 73 | public: | ||
| 74 | // We use convenient declarations that we derive from the property system | ||
| 75 | using CO2 = CO2Impl; | ||
| 76 | using H2O = H2OType; | ||
| 77 | // export the underlying brine fluid system for the liquid phase, as brine is used as a "pseudo component" | ||
| 78 | using Brine = Dumux::FluidSystems::ICPComplexSalinityBrine<Scalar, H2OType>; | ||
| 79 | using Na = Components::SodiumIon<Scalar>; | ||
| 80 | using Cl = Components::ChlorideIon<Scalar>; | ||
| 81 | using Ca = Components::CalciumIon<Scalar>; | ||
| 82 | using Urea = Components::Urea<Scalar>; | ||
| 83 | using O2 = Components::O2<Scalar>; | ||
| 84 | using Glucose = Components::Glucose<Scalar>; | ||
| 85 | using SuspendedBiomass = Components::SuspendedBiomass<Scalar>; | ||
| 86 | |||
| 87 | // We define the binary coefficients file, which accounts for the interactions of the main fluids in our setup, water/brine and CO2 | ||
| 88 | using Brine_CO2 = BinaryCoeff::Brine_CO2<Scalar, CO2Impl, true>; | ||
| 89 | |||
| 90 | // the type of parameter cache objects. this fluid system does not | ||
| 91 | // cache anything, so it uses Dumux::NullParameterCache | ||
| 92 | typedef Dumux::NullParameterCache ParameterCache; | ||
| 93 | |||
| 94 | // We define phase-related indices and properties | ||
| 95 | static constexpr int numPhases = 2; // liquid and gas phases | ||
| 96 | static constexpr int wPhaseIdx = 0; // index of the liquid phase | ||
| 97 | static constexpr int nPhaseIdx = 1; // index of the gas phase | ||
| 98 | static constexpr int phase0Idx = wPhaseIdx; | ||
| 99 | static constexpr int phase1Idx = nPhaseIdx; | ||
| 100 | |||
| 101 | // the phase names | ||
| 102 | 28 | static std::string phaseName(int phaseIdx) | |
| 103 | { | ||
| 104 |
8/16✓ Branch 0 taken 1 times.
✓ Branch 1 taken 27 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
28 | static std::string name[] = { |
| 105 | "w", | ||
| 106 | "n" | ||
| 107 | }; | ||
| 108 | |||
| 109 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 28 times.
|
28 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
| 110 | 28 | return name[phaseIdx]; | |
| 111 | } | ||
| 112 | |||
| 113 | // We define component-related indices and properties | ||
| 114 | static constexpr int numComponents = 9; // H2O, TotalC, Na, Cl, Ca,... | ||
| 115 | static constexpr int numSecComponents = 6; | ||
| 116 | |||
| 117 | static constexpr int H2OIdx = 0; | ||
| 118 | static constexpr int BrineIdx = 0; | ||
| 119 | static constexpr int TCIdx = 1; | ||
| 120 | static constexpr int wCompIdx = BrineIdx; | ||
| 121 | static constexpr int nCompIdx = TCIdx; | ||
| 122 | static constexpr int comp0Idx = wCompIdx; | ||
| 123 | static constexpr int comp1Idx = nCompIdx; | ||
| 124 | |||
| 125 | static constexpr int NaIdx = 2; | ||
| 126 | static constexpr int ClIdx = 3; | ||
| 127 | static constexpr int CaIdx = 4; | ||
| 128 | static constexpr int UreaIdx = 5; | ||
| 129 | static constexpr int O2Idx = 6; | ||
| 130 | static constexpr int GlucoseIdx = 7; | ||
| 131 | static constexpr int SuspendedBiomassIdx = 8; | ||
| 132 | // [[/codeblock]] | ||
| 133 | // [[exclude]] | ||
| 134 | |||
| 135 | // We check whether a phase is liquid | ||
| 136 | static constexpr bool isLiquid(int phaseIdx) | ||
| 137 | { | ||
| 138 | assert(0 <= phaseIdx && phaseIdx < numPhases); | ||
| 139 | |||
| 140 | return phaseIdx != nPhaseIdx; | ||
| 141 | } | ||
| 142 | |||
| 143 | // We return whether a phase is gaseous | ||
| 144 | static constexpr bool isGas(int phaseIdx) | ||
| 145 | { | ||
| 146 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5191544 times.
|
5191544 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
| 147 | return phaseIdx == nPhaseIdx; | ||
| 148 | } | ||
| 149 | |||
| 150 | // We assume ideal mixtures | ||
| 151 | static constexpr bool isIdealMixture(int phaseIdx) | ||
| 152 | { | ||
| 153 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 2595772 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2595772 times.
|
5191544 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
| 154 | |||
| 155 | return true; | ||
| 156 | } | ||
| 157 | |||
| 158 | // We assume compressible fluid phases | ||
| 159 | static constexpr bool isCompressible(int phaseIdx) | ||
| 160 | { | ||
| 161 | assert(0 <= phaseIdx && phaseIdx < numPhases); | ||
| 162 | |||
| 163 | return true; | ||
| 164 | } | ||
| 165 | |||
| 166 | // The component names | ||
| 167 | 18 | static std::string componentName(int compIdx) | |
| 168 | { | ||
| 169 | |||
| 170 |
9/10✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
|
18 | switch (compIdx) { |
| 171 | 2 | case BrineIdx: return H2O::name(); | |
| 172 | 4 | case TCIdx: return "TotalC"; | |
| 173 | 2 | case CaIdx: return Ca::name(); | |
| 174 | 2 | case NaIdx: return Na::name(); | |
| 175 | 2 | case ClIdx: return Cl::name(); | |
| 176 | 2 | case SuspendedBiomassIdx: return SuspendedBiomass::name(); | |
| 177 | 2 | case GlucoseIdx: return Glucose::name(); | |
| 178 | 2 | case O2Idx: return O2::name(); | |
| 179 | 2 | case UreaIdx: return Urea::name(); | |
| 180 | ✗ | default: DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); break; | |
| 181 | }; | ||
| 182 | } | ||
| 183 | |||
| 184 | // The component molar masses | ||
| 185 | 737396367 | static Scalar molarMass(int compIdx) | |
| 186 | { | ||
| 187 |
9/10✓ Branch 0 taken 54853913 times.
✓ Branch 1 taken 113415126 times.
✓ Branch 2 taken 113435711 times.
✓ Branch 3 taken 113414436 times.
✓ Branch 4 taken 67023443 times.
✓ Branch 5 taken 67032896 times.
✓ Branch 6 taken 67032896 times.
✓ Branch 7 taken 58912309 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 82275637 times.
|
737396367 | switch (compIdx) { |
| 188 | case H2OIdx: return H2O::molarMass(); | ||
| 189 | // actually, the molar mass of brine is only needed for diffusion | ||
| 190 | // but since chloride and sodium are accounted for separately | ||
| 191 | // only the molar mass of water is returned. | ||
| 192 | 54853913 | case TCIdx: return CO2::molarMass(); | |
| 193 | 113415126 | case CaIdx: return Ca::molarMass(); | |
| 194 | 113435711 | case NaIdx: return Na::molarMass(); | |
| 195 | 113414436 | case ClIdx: return Cl::molarMass(); | |
| 196 | 67023443 | case SuspendedBiomassIdx: return SuspendedBiomass::molarMass(); | |
| 197 | 67032896 | case GlucoseIdx: return Glucose::molarMass(); | |
| 198 | 67032896 | case O2Idx: return O2::molarMass(); | |
| 199 | 58912309 | case UreaIdx: return Urea::molarMass(); | |
| 200 | ✗ | default: DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); break; | |
| 201 | }; | ||
| 202 | } | ||
| 203 | // [[/exclude]] | ||
| 204 | // | ||
| 205 | // #### The Brine Adapter | ||
| 206 | // With the brine adapter, we link water and the components sodium, chloride, and calcium to form brine, to be able to use the "brine.hh" fluidsystem expecting only water and a single salt. | ||
| 207 | // Here, we define that the components water, sodium, chloride, and calcium contribute to brine. | ||
| 208 | // The real work is done by the adapter "icpcomplexsalinitybrine.hh". | ||
| 209 | // [[codeblock]] | ||
| 210 | private: | ||
| 211 | struct BrineAdapterPolicy | ||
| 212 | { | ||
| 213 | using FluidSystem = Brine; | ||
| 214 | |||
| 215 | ✗ | static constexpr int phaseIdx(int brinePhaseIdx) { return wPhaseIdx; } | |
| 216 | 137575916 | static constexpr int compIdx(int brineCompIdx) | |
| 217 | { | ||
| 218 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 137575916 times.
|
137575916 | assert( brineCompIdx == Brine::H2OIdx || brineCompIdx == Brine::NaIdx |
| 219 | || brineCompIdx == Brine::ClIdx || brineCompIdx == Brine::CaIdx); | ||
| 220 |
2/2✓ Branch 0 taken 28553492 times.
✓ Branch 1 taken 109022424 times.
|
137575916 | switch (brineCompIdx) |
| 221 | { | ||
| 222 | case Brine::H2OIdx: return H2OIdx; | ||
| 223 | case Brine::NaIdx: return NaIdx; | ||
| 224 | case Brine::ClIdx: return ClIdx; | ||
| 225 | case Brine::CaIdx: return CaIdx; | ||
| 226 | default: return 0; // this will never be reached, only needed to suppress compiler warning | ||
| 227 | } | ||
| 228 | } | ||
| 229 | }; | ||
| 230 | |||
| 231 | template<class FluidState> | ||
| 232 | using BrineAdapter = FluidStateAdapter<FluidState, BrineAdapterPolicy>; | ||
| 233 | // [[/codeblock]] | ||
| 234 | |||
| 235 | // #### Initializing the fluid system | ||
| 236 | // [[codeblock]] | ||
| 237 | public: | ||
| 238 | static void init() | ||
| 239 | { | ||
| 240 | init(/*startTemp=*/275.15, /*endTemp=*/455.15, /*tempSteps=*/30, | ||
| 241 | /*startPressure=*/1e4, /*endPressure=*/99e6, /*pressureSteps=*/500); | ||
| 242 | } | ||
| 243 | static void init(Scalar startTemp, Scalar endTemp, int tempSteps, | ||
| 244 | Scalar startPressure, Scalar endPressure, int pressureSteps) | ||
| 245 | { | ||
| 246 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | std::cout << "Initializing tables for the pure-water properties.\n"; |
| 247 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | H2O::init(startTemp, endTemp, tempSteps, |
| 248 | startPressure, endPressure, pressureSteps); | ||
| 249 | } | ||
| 250 | // [[/codeblock]] | ||
| 251 | |||
| 252 | // #### The Component-Phase Interactions | ||
| 253 | // In the following, we define the component-phase interactions such as // each component's fugacity coefficient in each phase | ||
| 254 | // and each component's and the phases' main components binary diffusion coefficient in the respective phase. | ||
| 255 | // [[codeblock]] | ||
| 256 | // The component fugacity coefficients | ||
| 257 | template <class FluidState> | ||
| 258 | 46723896 | static Scalar fugacityCoefficient(const FluidState &fluidState, | |
| 259 | const ParameterCache ¶mCache, | ||
| 260 | int phaseIdx, | ||
| 261 | int compIdx) | ||
| 262 | { | ||
| 263 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 46723896 times.
|
46723896 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
| 264 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 46723896 times.
|
46723896 | assert(0 <= compIdx && compIdx < numComponents); |
| 265 | |||
| 266 |
2/2✓ Branch 0 taken 23361948 times.
✓ Branch 1 taken 23361948 times.
|
46723896 | if (phaseIdx == nPhaseIdx) |
| 267 | // use the fugacity coefficients of an ideal gas. the | ||
| 268 | // actual value of the fugacity is not relevant, as long | ||
| 269 | // as the relative fluid compositions are observed, | ||
| 270 | return 1.0; | ||
| 271 | |||
| 272 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23361948 times.
|
23361948 | Scalar temperature = fluidState.temperature(phaseIdx); |
| 273 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23361948 times.
|
23361948 | Scalar pressure = fluidState.pressure(phaseIdx); |
| 274 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23361948 times.
|
23361948 | if (pressure<0) |
| 275 | { | ||
| 276 | typedef Dune::FieldVector<Scalar, numComponents> ComponentVector; | ||
| 277 | ✗ | ComponentVector moleFractionsw; | |
| 278 | ✗ | ComponentVector massFractionsw; | |
| 279 | |||
| 280 | ✗ | for (int compIdx = 0; compIdx<numComponents;++compIdx) | |
| 281 | { | ||
| 282 | ✗ | moleFractionsw[compIdx] = fluidState.moleFraction(wPhaseIdx,compIdx); | |
| 283 | ✗ | massFractionsw[compIdx] = fluidState.massFraction(wPhaseIdx,compIdx); | |
| 284 | } | ||
| 285 | } | ||
| 286 | |||
| 287 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23361948 times.
|
23361948 | assert(temperature > 0); |
| 288 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23361948 times.
|
23361948 | assert(pressure > 0); |
| 289 | |||
| 290 | // calculate the equilibrium composition for the given | ||
| 291 | // temperature and pressure. | ||
| 292 | Scalar xgH2O, xlH2O; | ||
| 293 | Scalar xlCO2, xgCO2; | ||
| 294 | 23361948 | Scalar Xl_Sal = fluidState.massFraction(wPhaseIdx, NaIdx) //Salinity= XNa+XCl+XCa | |
| 295 | 23361948 | + fluidState.massFraction(wPhaseIdx, ClIdx) | |
| 296 | 23361948 | + fluidState.massFraction(wPhaseIdx, CaIdx); | |
| 297 | 23361948 | Brine_CO2::calculateMoleFractions(temperature, | |
| 298 | pressure, | ||
| 299 | Xl_Sal, | ||
| 300 | /*knownPhaseIdx=*/ -1, | ||
| 301 | xlCO2, | ||
| 302 | xgH2O); | ||
| 303 | |||
| 304 |
2/4✓ Branch 0 taken 23361948 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 23361948 times.
✗ Branch 3 not taken.
|
46723896 | xlCO2 = std::max(0.0, std::min(1.0, xlCO2)); |
| 305 |
2/4✓ Branch 0 taken 23361948 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 23361948 times.
✗ Branch 3 not taken.
|
46723896 | xgH2O = std::max(0.0, std::min(1.0, xgH2O)); |
| 306 | 23361948 | xlH2O = 1.0 - xlCO2; | |
| 307 | 23361948 | xgCO2 = 1.0 - xgH2O; | |
| 308 | |||
| 309 | // Only H2O, CO2, and O2 in the gas phase | ||
| 310 |
2/2✓ Branch 0 taken 2595772 times.
✓ Branch 1 taken 20766176 times.
|
23361948 | if (compIdx == BrineIdx) { |
| 311 | 2595772 | Scalar phigH2O = 1.0; | |
| 312 | 2595772 | return phigH2O * xgH2O / xlH2O; | |
| 313 | } | ||
| 314 |
2/2✓ Branch 0 taken 2595772 times.
✓ Branch 1 taken 18170404 times.
|
20766176 | else if (compIdx == TCIdx) |
| 315 | { | ||
| 316 | 2595772 | Scalar phigCO2 = 1.0; | |
| 317 | 2595772 | return phigCO2 * xgCO2 / xlCO2; | |
| 318 | } | ||
| 319 |
2/2✓ Branch 0 taken 2595772 times.
✓ Branch 1 taken 15574632 times.
|
18170404 | else if (compIdx == O2Idx) |
| 320 | { | ||
| 321 | 2595772 | return Dumux::BinaryCoeff::H2O_O2::henry(temperature)/pressure; | |
| 322 | } | ||
| 323 | // All other components are assumed not be present in the gas phase | ||
| 324 | else | ||
| 325 | { | ||
| 326 | return 1e-20; | ||
| 327 | } | ||
| 328 | } | ||
| 329 | |||
| 330 | template <class FluidState> | ||
| 331 | static Scalar diffusionCoefficient(const FluidState &fluidState, | ||
| 332 | const ParameterCache ¶mCache, | ||
| 333 | int phaseIdx, | ||
| 334 | int compIdx) | ||
| 335 | { | ||
| 336 | DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients"); | ||
| 337 | } | ||
| 338 | |||
| 339 | // The binary diffusion coefficients of the components and the phases main component | ||
| 340 | template <class FluidState> | ||
| 341 | 41532352 | static Scalar binaryDiffusionCoefficient(const FluidState &fluidState, | |
| 342 | const ParameterCache ¶mCache, | ||
| 343 | int phaseIdx, | ||
| 344 | int compIIdx, | ||
| 345 | int compJIdx) | ||
| 346 | { | ||
| 347 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 41532352 times.
|
41532352 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
| 348 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 41532352 times.
|
41532352 | assert(0 <= compIIdx && compIIdx < numComponents); |
| 349 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 41532352 times.
|
41532352 | assert(0 <= compJIdx && compJIdx < numComponents); |
| 350 | |||
| 351 |
2/2✓ Branch 0 taken 20766176 times.
✓ Branch 1 taken 20766176 times.
|
41532352 | const Scalar temperature = fluidState.temperature(phaseIdx); |
| 352 |
2/2✓ Branch 0 taken 20766176 times.
✓ Branch 1 taken 20766176 times.
|
41532352 | const Scalar pressure = fluidState.pressure(phaseIdx); |
| 353 | |||
| 354 |
2/2✓ Branch 0 taken 20766176 times.
✓ Branch 1 taken 20766176 times.
|
41532352 | if (phaseIdx == wPhaseIdx) |
| 355 | { | ||
| 356 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 20766176 times.
|
20766176 | assert(compIIdx == H2OIdx); |
| 357 |
2/2✓ Branch 0 taken 2595772 times.
✓ Branch 1 taken 18170404 times.
|
20766176 | if (compJIdx == TCIdx) |
| 358 | 2595772 | return Brine_CO2::liquidDiffCoeff(temperature, pressure); | |
| 359 |
2/2✓ Branch 0 taken 2595772 times.
✓ Branch 1 taken 15574632 times.
|
18170404 | else if(compJIdx == O2Idx) |
| 360 | 5191544 | return Dumux::BinaryCoeff::H2O_O2::liquidDiffCoeff(temperature, pressure); | |
| 361 | else //all other components | ||
| 362 | return 1.587e-9; //[m²/s] //J. Phys. D: Appl. Phys. 40 (2007) 2769-2776 //old Value from Anozie 1e-9 | ||
| 363 | } | ||
| 364 | else | ||
| 365 | { | ||
| 366 | assert(phaseIdx == nPhaseIdx); | ||
| 367 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 20766176 times.
|
20766176 | assert(compIIdx == TCIdx); |
| 368 |
2/2✓ Branch 0 taken 2595772 times.
✓ Branch 1 taken 18170404 times.
|
20766176 | if (compJIdx == H2OIdx) |
| 369 | 2595772 | return Brine_CO2::gasDiffCoeff(temperature, pressure); | |
| 370 |
2/2✓ Branch 0 taken 2595772 times.
✓ Branch 1 taken 15574632 times.
|
18170404 | else if (compJIdx == O2Idx) |
| 371 | 2595772 | return Dumux::BinaryCoeff::H2O_O2::gasDiffCoeff(temperature, pressure); | |
| 372 | else //all other components | ||
| 373 | return 0.0; | ||
| 374 | } | ||
| 375 | } | ||
| 376 | // [[/codeblock]] | ||
| 377 | |||
| 378 | // #### The Fluid Properties | ||
| 379 | // In the following, all functions defining the phase properties are given: | ||
| 380 | // the density, viscosity, enthalpy, thermal conductivities, and heat capacities | ||
| 381 | // of each phase depending on temperature, pressure, and phase composition | ||
| 382 | // [[codeblock]] | ||
| 383 | // The phase density of the liquid phase is calculated accounting for the impact of solutes in the brine as well as the contribution of CO2. | ||
| 384 | // For the gas phase, a mixture of water and CO2 is considered, as solutes do not partition into the gas phase. | ||
| 385 | template <class FluidState> | ||
| 386 | 5191544 | static Scalar density(const FluidState &fluidState, | |
| 387 | const ParameterCache ¶mCache, | ||
| 388 | int phaseIdx) | ||
| 389 | { | ||
| 390 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5191544 times.
|
5191544 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
| 391 | |||
| 392 |
2/2✓ Branch 0 taken 2595772 times.
✓ Branch 1 taken 2595772 times.
|
5191544 | const Scalar T = fluidState.temperature(phaseIdx); |
| 393 | |||
| 394 |
2/2✓ Branch 0 taken 2595772 times.
✓ Branch 1 taken 2595772 times.
|
5191544 | if (phaseIdx == wPhaseIdx) |
| 395 | { | ||
| 396 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2595772 times.
|
2595772 | const Scalar pl = fluidState.pressure(phaseIdx); |
| 397 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2595772 times.
|
2595772 | const Scalar xlCO2 = fluidState.moleFraction(wPhaseIdx, TCIdx); |
| 398 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2595772 times.
|
2595772 | const Scalar xlH2O = fluidState.moleFraction(wPhaseIdx, H2OIdx); |
| 399 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2595772 times.
|
2595772 | if(T < 273.15) { |
| 400 | ✗ | DUNE_THROW(NumericalProblem, | |
| 401 | "Liquid density for Brine and CO2 is only " | ||
| 402 | "defined above 273.15K (is" << T << ")"); | ||
| 403 | } | ||
| 404 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2595772 times.
|
2595772 | if(pl >= 2.5e8) { |
| 405 | ✗ | DUNE_THROW(NumericalProblem, | |
| 406 | "Liquid density for Brine and CO2 is only " | ||
| 407 | "defined below 250MPa (is" << pl << ")"); | ||
| 408 | } | ||
| 409 | |||
| 410 | // calling the brine adapter for brine density | ||
| 411 | 5191544 | Scalar rho_brine = Brine::molarDensity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx) | |
| 412 | 2595772 | *(H2O::molarMass()*fluidState.moleFraction(wPhaseIdx, H2OIdx) | |
| 413 | 2595772 | + Na::molarMass()*fluidState.moleFraction(wPhaseIdx, NaIdx) | |
| 414 | 2595772 | + Cl::molarMass()*fluidState.moleFraction(wPhaseIdx, ClIdx) | |
| 415 | 2595772 | + Ca::molarMass()*fluidState.moleFraction(wPhaseIdx, CaIdx) | |
| 416 | ); | ||
| 417 | // the density of pure water | ||
| 418 | 2595772 | Scalar rho_pure = H2O::liquidDensity(T, pl); | |
| 419 | //we also use a private helper function to calculate the liquid density of the same amount of CO2 in pure water | ||
| 420 | 2595772 | Scalar rho_lCO2 = liquidDensityWaterCO2_(T, pl, xlH2O, xlCO2); | |
| 421 | // calculating the density contribution of CO2 in pure water | ||
| 422 | 2595772 | Scalar contribCO2 = rho_lCO2 - rho_pure; | |
| 423 | // assuming CO2 increases the density for brine by the same amount as for pure water | ||
| 424 | 2595772 | return rho_brine + contribCO2; | |
| 425 | } | ||
| 426 | |||
| 427 | else if (phaseIdx == nPhaseIdx) | ||
| 428 | 5191544 | return H2O::gasDensity(T, fluidState.partialPressure(nPhaseIdx, H2OIdx)) | |
| 429 | 7787316 | + CO2::gasDensity(T, fluidState.partialPressure(nPhaseIdx, TCIdx)); | |
| 430 | else | ||
| 431 | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | ||
| 432 | } | ||
| 433 | |||
| 434 | // The phase molar density of the liquid phase is assumed to not change significantly from the molar density of the pure brine. | ||
| 435 | // For the gas phase, a mixture of water and CO2 is considered, as solutes do not partition into the gas phase. | ||
| 436 | using Base::molarDensity; | ||
| 437 | template <class FluidState> | ||
| 438 | 5191544 | static Scalar molarDensity(const FluidState& fluidState, int phaseIdx) | |
| 439 | { | ||
| 440 |
2/2✓ Branch 0 taken 2595772 times.
✓ Branch 1 taken 2595772 times.
|
5191544 | if (phaseIdx == wPhaseIdx) |
| 441 | 7787316 | return Brine::molarDensity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx); | |
| 442 |
1/2✓ Branch 0 taken 2595772 times.
✗ Branch 1 not taken.
|
2595772 | else if (phaseIdx == nPhaseIdx) |
| 443 | 2595772 | return H2O::gasMolarDensity(fluidState.temperature(phaseIdx), | |
| 444 | fluidState.partialPressure(phaseIdx, H2OIdx)) | ||
| 445 | 2595772 | + CO2::gasMolarDensity(fluidState.temperature(phaseIdx), | |
| 446 | 2595772 | fluidState.partialPressure(phaseIdx, TCIdx)); | |
| 447 | else | ||
| 448 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | |
| 449 | } | ||
| 450 | |||
| 451 | // The phase viscosities are assumed to not change significantly from those of the pure brine for the liquid and pure CO2 for the gas phase | ||
| 452 | using Base::viscosity; | ||
| 453 | template <class FluidState> | ||
| 454 | 5191544 | static Scalar viscosity(const FluidState& fluidState, int phaseIdx) | |
| 455 | { | ||
| 456 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5191544 times.
|
5191544 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
| 457 | |||
| 458 |
2/2✓ Branch 0 taken 2595772 times.
✓ Branch 1 taken 2595772 times.
|
5191544 | if (phaseIdx == wPhaseIdx) |
| 459 | 5191544 | return Brine::viscosity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx); | |
| 460 | else if (phaseIdx == nPhaseIdx) | ||
| 461 | { | ||
| 462 | 7787316 | return CO2::gasViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); | |
| 463 | } | ||
| 464 | else | ||
| 465 | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | ||
| 466 | |||
| 467 | } | ||
| 468 | // [[/codeblock]] | ||
| 469 | |||
| 470 | // [[exclude]] | ||
| 471 | // non-isothermal phase properties | ||
| 472 | // The phase enthalpies are assumed to not change significantly from those of the pure brine for the liquid and pure CO2 for the gas phase | ||
| 473 | template <class FluidState> | ||
| 474 | static Scalar enthalpy(const FluidState &fluidState, | ||
| 475 | const ParameterCache ¶mCache, | ||
| 476 | int phaseIdx) | ||
| 477 | { | ||
| 478 | assert(0 <= phaseIdx && phaseIdx < numPhases); | ||
| 479 | |||
| 480 | const Scalar temperature = fluidState.temperature(phaseIdx); | ||
| 481 | const Scalar pressure = fluidState.pressure(phaseIdx); | ||
| 482 | |||
| 483 | if (phaseIdx == wPhaseIdx) | ||
| 484 | return Brine::enthalpy(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx); | ||
| 485 | else | ||
| 486 | return fluidState.massFraction(nPhaseIdx, H2OIdx)*Brine::gasEnthalpy(temperature, pressure) | ||
| 487 | + fluidState.massFraction(nPhaseIdx, TCIdx)*CO2::gasEnthalpy(temperature, pressure); | ||
| 488 | } | ||
| 489 | |||
| 490 | // The phase thermal conductivities are assumed to not change significantly from those of the pure brine for the liquid and pure CO2 for the gas phase | ||
| 491 | template <class FluidState> | ||
| 492 | static Scalar thermalConductivity(const FluidState& fluidState, int phaseIdx) | ||
| 493 | { | ||
| 494 | if (phaseIdx == wPhaseIdx) | ||
| 495 | return Brine::thermalConductivity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx); | ||
| 496 | else if (phaseIdx ==nPhaseIdx) | ||
| 497 | return CO2::gasThermalConductivity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); | ||
| 498 | |||
| 499 | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | ||
| 500 | } | ||
| 501 | |||
| 502 | // The phase heat capacitiy of the liquid phase is assumed to not change significantly from the heat capacitiy of the pure brine. | ||
| 503 | // For the gas phase, a mixture of water and CO2 is considered, as solutes do not partition into the gas phase. | ||
| 504 | template <class FluidState> | ||
| 505 | static Scalar heatCapacity(const FluidState &fluidState, int phaseIdx) | ||
| 506 | { | ||
| 507 | const Scalar T = fluidState.temperature(phaseIdx); | ||
| 508 | const Scalar p = fluidState.pressure(phaseIdx); | ||
| 509 | |||
| 510 | if (phaseIdx == wPhaseIdx) | ||
| 511 | return Brine::heatCapacity(BrineAdapter<FluidState>(fluidState), Brine::liquidPhaseIdx); | ||
| 512 | |||
| 513 | // We assume NaCl not to be present in the gas phase here | ||
| 514 | else if (phaseIdx ==nPhaseIdx) | ||
| 515 | return CO2::gasHeatCapacity(T, p)*fluidState.moleFraction(nPhaseIdx, TCIdx) | ||
| 516 | + H2O::gasHeatCapacity(T, p)*fluidState.moleFraction(nPhaseIdx, H2OIdx); | ||
| 517 | |||
| 518 | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | ||
| 519 | } | ||
| 520 | |||
| 521 | |||
| 522 | private: | ||
| 523 | |||
| 524 | //a private helper function to calculate the liquid density of the same amount of CO2 as we have in our brine in pure water | ||
| 525 | ✗ | static Scalar liquidDensityWaterCO2_(Scalar temperature, | |
| 526 | Scalar pl, | ||
| 527 | Scalar xlH2O, | ||
| 528 | Scalar xlCO2) | ||
| 529 | { | ||
| 530 | ✗ | const Scalar M_CO2 = CO2::molarMass(); | |
| 531 | ✗ | const Scalar M_H2O = H2O::molarMass(); | |
| 532 | |||
| 533 | ✗ | const Scalar tempC = temperature - 273.15; /* tempC : temperature in °C */ | |
| 534 | ✗ | const Scalar rho_pure = H2O::liquidDensity(temperature, pl); | |
| 535 | ✗ | xlH2O = 1.0 - xlCO2; // xlH2O is available, but in case of a pure gas phase | |
| 536 | // the value of M_T for the virtual liquid phase can become very large | ||
| 537 | ✗ | const Scalar M_T = M_H2O * xlH2O + M_CO2 * xlCO2; | |
| 538 | ✗ | const Scalar V_phi = | |
| 539 | ✗ | (37.51 + | |
| 540 | ✗ | tempC*(-9.585e-2 + | |
| 541 | ✗ | tempC*(8.74e-4 - | |
| 542 | ✗ | tempC*5.044e-7))) / 1.0e6; | |
| 543 | ✗ | return 1/ (xlCO2 * V_phi/M_T + M_H2O * xlH2O / (rho_pure * M_T)); | |
| 544 | } | ||
| 545 | }; | ||
| 546 | |||
| 547 | } // end namespace Dumux::FluidSystems | ||
| 548 | // [[/exclude]] | ||
| 549 | // [[/content]] | ||
| 550 | #endif | ||
| 551 |