3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
porousmediumflow/2pnc/volumevariables.hh
Go to the documentation of this file.
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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
26#ifndef DUMUX_2PNC_VOLUME_VARIABLES_HH
27#define DUMUX_2PNC_VOLUME_VARIABLES_HH
28
29#include <iostream>
30#include <vector>
31
32#include <dumux/common/math.hh>
34
37
42
44
46
47namespace Dumux {
48
54template <class Traits>
57, public EnergyVolumeVariables<Traits, TwoPNCVolumeVariables<Traits> >
58{
61 using Scalar = typename Traits::PrimaryVariables::value_type;
62 using PermeabilityType = typename Traits::PermeabilityType;
63 using FS = typename Traits::FluidSystem;
64 using ModelTraits = typename Traits::ModelTraits;
65 static constexpr int numFluidComps = ParentType::numFluidComponents();
66 enum
67 {
68 numMajorComponents = ModelTraits::numFluidPhases(),
69
70 // phase indices
71 phase0Idx = FS::phase0Idx,
72 phase1Idx = FS::phase1Idx,
73
74 // component indices
75 comp0Idx = FS::comp0Idx,
76 comp1Idx = FS::comp1Idx,
77
78 // phase presence enums
79 secondPhaseOnly = ModelTraits::Indices::secondPhaseOnly,
80 firstPhaseOnly = ModelTraits::Indices::firstPhaseOnly,
81 bothPhases = ModelTraits::Indices::bothPhases,
82
83 // primary variable indices
84 pressureIdx = ModelTraits::Indices::pressureIdx,
85 switchIdx = ModelTraits::Indices::switchIdx
86 };
87
88 static constexpr auto formulation = ModelTraits::priVarFormulation();
89 static constexpr bool setFirstPhaseMoleFractions = ModelTraits::setMoleFractionsForFirstPhase();
90
93 using EffDiffModel = typename Traits::EffectiveDiffusivityModel;
94 using DiffusionCoefficients = typename Traits::DiffusionType::DiffusionCoefficientsContainer;
95
96public:
98 using FluidState = typename Traits::FluidState;
100 using FluidSystem = typename Traits::FluidSystem;
102 using Indices = typename ModelTraits::Indices;
104 using SolidState = typename Traits::SolidState;
106 using SolidSystem = typename Traits::SolidSystem;
109
111 static constexpr bool useMoles() { return Traits::ModelTraits::useMoles(); }
113 static constexpr TwoPFormulation priVarFormulation() { return formulation; }
114
115 // check for permissive specifications
116 static_assert(useMoles(), "use moles has to be set true in the 2pnc model");
117 static_assert(ModelTraits::numFluidPhases() == 2, "NumPhases set in the model is not two!");
118 static_assert((formulation == TwoPFormulation::p0s1 || formulation == TwoPFormulation::p1s0), "Chosen TwoPFormulation not supported!");
119
129 template<class ElemSol, class Problem, class Element, class Scv>
130 void update(const ElemSol &elemSol,
131 const Problem &problem,
132 const Element &element,
133 const Scv& scv)
134 {
135 ParentType::update(elemSol, problem, element, scv);
136
137 completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
138
140 // calculate the remaining quantities
142
143 const auto& spatialParams = problem.spatialParams();
144 const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol);
145
146 const int wPhaseIdx = fluidState_.wettingPhase();
147 const int nPhaseIdx = 1 - wPhaseIdx;
148
149 // mobilities -> require wetting phase saturation as parameter!
150 mobility_[wPhaseIdx] = fluidMatrixInteraction.krw(saturation(wPhaseIdx))/fluidState_.viscosity(wPhaseIdx);
151 mobility_[nPhaseIdx] = fluidMatrixInteraction.krn(saturation(wPhaseIdx))/fluidState_.viscosity(nPhaseIdx);
152
153 //update porosity before calculating the effective properties depending on it
154 updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
155
156 auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
157 {
158 return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
159 };
160
161 effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
162
163 // calculate the remaining quantities
164 EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
165 permeability_ = spatialParams.permeability(element, scv, elemSol);
166 EnergyVolVars::updateEffectiveThermalConductivity();
167 }
168
181 template<class ElemSol, class Problem, class Element, class Scv>
182 void completeFluidState(const ElemSol& elemSol,
183 const Problem& problem,
184 const Element& element,
185 const Scv& scv,
188 {
189 EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState);
190
191 const auto& priVars = elemSol[scv.localDofIndex()];
192 const auto phasePresence = priVars.state();
193
194 const auto& spatialParams = problem.spatialParams();
195 const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol);
196 const auto wPhaseIdx = spatialParams.template wettingPhase<FluidSystem>(element, scv, elemSol);
197 fluidState.setWettingPhase(wPhaseIdx);
198
199 // set the saturations
200 if (phasePresence == secondPhaseOnly)
201 {
202 fluidState.setSaturation(phase0Idx, 0.0);
203 fluidState.setSaturation(phase1Idx, 1.0);
204 }
205 else if (phasePresence == firstPhaseOnly)
206 {
207 fluidState.setSaturation(phase0Idx, 1.0);
208 fluidState.setSaturation(phase1Idx, 0.0);
209 }
210 else if (phasePresence == bothPhases)
211 {
212 if (formulation == TwoPFormulation::p0s1)
213 {
214 fluidState.setSaturation(phase1Idx, priVars[switchIdx]);
215 fluidState.setSaturation(phase0Idx, 1 - priVars[switchIdx]);
216 }
217 else
218 {
219 fluidState.setSaturation(phase0Idx, priVars[switchIdx]);
220 fluidState.setSaturation(phase1Idx, 1 - priVars[switchIdx]);
221 }
222 }
223 else
224 DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
225
226 // set pressures of the fluid phases
227 pc_ = fluidMatrixInteraction.pc(fluidState.saturation(wPhaseIdx));
228 if (formulation == TwoPFormulation::p0s1)
229 {
230 fluidState.setPressure(phase0Idx, priVars[pressureIdx]);
231 fluidState.setPressure(phase1Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] + pc_
232 : priVars[pressureIdx] - pc_);
233 }
234 else
235 {
236 fluidState.setPressure(phase1Idx, priVars[pressureIdx]);
237 fluidState.setPressure(phase0Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] - pc_
238 : priVars[pressureIdx] + pc_);
239 }
240
241 // calculate the phase compositions
242 typename FluidSystem::ParameterCache paramCache;
243
244 // now comes the tricky part: calculate phase composition
245 if (phasePresence == bothPhases)
246 {
247 // both phases are present, phase composition results from
248 // the first <-> second phase equilibrium. This is the job
249 // of the "MiscibleMultiPhaseComposition" constraint solver
250
251 // set the known mole fractions in the fluidState so that they
252 // can be used by the MiscibleMultiPhaseComposition constraint solver
253
254 const int knownPhaseIdx = setFirstPhaseMoleFractions ? phase0Idx : phase1Idx;
255 for (int compIdx = numMajorComponents; compIdx < ModelTraits::numFluidComponents(); ++compIdx)
256 fluidState.setMoleFraction(knownPhaseIdx, compIdx, priVars[compIdx]);
257
259 paramCache,
260 knownPhaseIdx);
261 }
262 else if (phasePresence == secondPhaseOnly)
263 {
264 Dune::FieldVector<Scalar, ModelTraits::numFluidComponents()> moleFrac;
265
266 moleFrac[comp0Idx] = priVars[switchIdx];
267 Scalar sumMoleFracOtherComponents = moleFrac[comp0Idx];
268
269 for (int compIdx = numMajorComponents; compIdx < ModelTraits::numFluidComponents(); ++compIdx)
270 {
271 moleFrac[compIdx] = priVars[compIdx];
272 sumMoleFracOtherComponents += moleFrac[compIdx];
273 }
274
275 moleFrac[comp1Idx] = 1 - sumMoleFracOtherComponents;
276
277 // Set fluid state mole fractions
278 for (int compIdx = 0; compIdx < ModelTraits::numFluidComponents(); ++compIdx)
279 fluidState.setMoleFraction(phase1Idx, compIdx, moleFrac[compIdx]);
280
281 // calculate the composition of the remaining phases (as
282 // well as the densities of all phases). this is the job
283 // of the "ComputeFromReferencePhase" constraint solver
285 paramCache,
286 phase1Idx);
287 }
288 else if (phasePresence == firstPhaseOnly)
289 {
290 // only the first phase is present, i.e. first phase composition
291 // is stored explicitly. extract _mass_ fractions in the second phase
292 Dune::FieldVector<Scalar, ModelTraits::numFluidComponents()> moleFrac;
293
294 moleFrac[comp1Idx] = priVars[switchIdx];
295 Scalar sumMoleFracOtherComponents = moleFrac[comp1Idx];
296 for (int compIdx = numMajorComponents; compIdx < ModelTraits::numFluidComponents(); ++compIdx)
297 {
298 moleFrac[compIdx] = priVars[compIdx];
299
300 sumMoleFracOtherComponents += moleFrac[compIdx];
301 }
302
303 moleFrac[comp0Idx] = 1 - sumMoleFracOtherComponents;
304
305 // convert mass to mole fractions and set the fluid state
306 for (int compIdx = 0; compIdx < ModelTraits::numFluidComponents(); ++compIdx)
307 fluidState.setMoleFraction(phase0Idx, compIdx, moleFrac[compIdx]);
308
309 // calculate the composition of the remaining phases (as
310 // well as the densities of all phases). this is the job
311 // of the "ComputeFromReferencePhase" constraint solver
313 paramCache,
314 phase0Idx);
315 }
316 paramCache.updateAll(fluidState);
317 for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
318 {
319 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
320 Scalar rhoMolar = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
321 Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
322 Scalar h = EnergyVolVars::enthalpy(fluidState, paramCache, phaseIdx);
323
324 fluidState.setDensity(phaseIdx, rho);
325 fluidState.setMolarDensity(phaseIdx, rhoMolar);
326 fluidState.setViscosity(phaseIdx, mu);
327 fluidState.setEnthalpy(phaseIdx, h);
328 }
329 }
330
334 const FluidState &fluidState() const
335 { return fluidState_; }
336
340 const SolidState &solidState() const
341 { return solidState_; }
342
348 Scalar averageMolarMass(int phaseIdx) const
349 { return fluidState_.averageMolarMass(phaseIdx); }
350
357 Scalar saturation(int phaseIdx) const
358 { return fluidState_.saturation(phaseIdx); }
359
366 Scalar density(int phaseIdx) const
367 { return fluidState_.density(phaseIdx); }
368
375 Scalar viscosity(int phaseIdx) const
376 { return fluidState_.viscosity(phaseIdx); }
377
384 Scalar molarDensity(int phaseIdx) const
385 {
386 if (phaseIdx < ModelTraits::numFluidPhases())
387 return fluidState_.molarDensity(phaseIdx);
388
389 else
390 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
391 }
392
399 Scalar pressure(int phaseIdx) const
400 { return fluidState_.pressure(phaseIdx); }
401
409 Scalar temperature() const
410 { return fluidState_.temperature(/*phaseIdx=*/0); }
411
418 Scalar mobility(int phaseIdx) const
419 { return mobility_[phaseIdx]; }
420
425 Scalar capillaryPressure() const
426 { return pc_; }
427
431 Scalar porosity() const
432 { return solidState_.porosity(); }
433
437 const PermeabilityType& permeability() const
438 { return permeability_; }
439
443 Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
444 {
445 typename FluidSystem::ParameterCache paramCache;
446 paramCache.updatePhase(fluidState_, phaseIdx);
447 return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx);
448 }
449
453 Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
454 { return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); }
455
462 Scalar massFraction(int phaseIdx, int compIdx) const
463 { return fluidState_.massFraction(phaseIdx, compIdx); }
464
471 Scalar moleFraction(int phaseIdx, int compIdx) const
472 { return fluidState_.moleFraction(phaseIdx, compIdx); }
473
477 int wettingPhase() const
478 { return fluidState_.wettingPhase(); }
479
480protected:
483
484private:
485
486 Scalar pc_; // The capillary pressure
487 Scalar porosity_; // Effective porosity within the control volume
488 PermeabilityType permeability_; // Effective permeability within the control volume
489 Scalar mobility_[ModelTraits::numFluidPhases()]; // Effective mobility within the control volume
490 DiffusionCoefficients effectiveDiffCoeff_;
491
492};
493
494} // end namespace Dumux
495
496#endif
Define some often used mathematical functions.
The available discretization methods in Dumux.
Computes all quantities of a generic fluid state if a reference phase has been specified.
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Update the solid volume fractions (inert and reacitve) and set them in the solidstate.
Defines an enumeration for the formulations accepted by the two-phase model.
TwoPFormulation
Enumerates the formulations which the two-phase model accepts.
Definition: formulation.hh:35
@ p1s0
first phase saturation and second phase pressure as primary variables
@ p0s1
first phase pressure and second phase saturation as primary variables
void updateSolidVolumeFractions(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv, SolidState &solidState, const int solidVolFracOffset)
update the solid volume fractions (inert and reacitve) and set them in the solidstate
Definition: updatesolidvolumefractions.hh:36
Definition: adapt.hh:29
std::string phasePresence() noexcept
I/O name of phase presence.
Definition: name.hh:147
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
std::string molarDensity(int phaseIdx) noexcept
I/O name of molar density for multiphase systems.
Definition: name.hh:83
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Computes all quantities of a generic fluid state if a reference phase has been specified.
Definition: computefromreferencephase.hh:76
static void solve(FluidState &fluidState, ParameterCache &paramCache, int refPhaseIdx)
Computes all quantities of a generic fluid state if a reference phase has been specified.
Definition: computefromreferencephase.hh:111
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: misciblemultiphasecomposition.hh:59
static void solve(FluidState &fluidState, ParameterCache &paramCache, int knownPhaseIdx=0)
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: misciblemultiphasecomposition.hh:81
The primary variable switch controlling the phase presence state variable.
Definition: 2pnc/primaryvariableswitch.hh:41
Contains the quantities which are are constant within a finite volume in the two-phase,...
Definition: porousmediumflow/2pnc/volumevariables.hh:58
Scalar porosity() const
Returns the average porosity within the control volume.
Definition: porousmediumflow/2pnc/volumevariables.hh:431
Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
Returns the binary diffusion coefficients for a phase in .
Definition: porousmediumflow/2pnc/volumevariables.hh:443
Scalar molarDensity(int phaseIdx) const
Returns the molar density of a given phase within the control volume.
Definition: porousmediumflow/2pnc/volumevariables.hh:384
static constexpr TwoPFormulation priVarFormulation()
Return the two-phase formulation used here.
Definition: porousmediumflow/2pnc/volumevariables.hh:113
Scalar viscosity(int phaseIdx) const
Returns the kinematic viscosity of a given phase within the control volume.
Definition: porousmediumflow/2pnc/volumevariables.hh:375
Scalar capillaryPressure() const
Returns the effective capillary pressure within the control volume in .
Definition: porousmediumflow/2pnc/volumevariables.hh:425
Scalar density(int phaseIdx) const
Returns the mass density of a given phase within the control volume.
Definition: porousmediumflow/2pnc/volumevariables.hh:366
void update(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv)
Updates all quantities for a given control volume.
Definition: porousmediumflow/2pnc/volumevariables.hh:130
const SolidState & solidState() const
Returns the phase state for the control-volume.
Definition: porousmediumflow/2pnc/volumevariables.hh:340
Scalar moleFraction(int phaseIdx, int compIdx) const
Returns the mole fraction of a component in the phase.
Definition: porousmediumflow/2pnc/volumevariables.hh:471
typename ModelTraits::Indices Indices
Export the indices.
Definition: porousmediumflow/2pnc/volumevariables.hh:102
typename Traits::SolidState SolidState
Export type of solid state.
Definition: porousmediumflow/2pnc/volumevariables.hh:104
Scalar saturation(int phaseIdx) const
Returns the saturation of a given phase within the control volume in .
Definition: porousmediumflow/2pnc/volumevariables.hh:357
typename Traits::SolidSystem SolidSystem
Export type of solid system.
Definition: porousmediumflow/2pnc/volumevariables.hh:106
const PermeabilityType & permeability() const
Returns the permeability within the control volume.
Definition: porousmediumflow/2pnc/volumevariables.hh:437
Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
Returns the effective diffusion coefficients for a phase in .
Definition: porousmediumflow/2pnc/volumevariables.hh:453
Scalar temperature() const
Returns temperature inside the sub-control volume.
Definition: porousmediumflow/2pnc/volumevariables.hh:409
Scalar averageMolarMass(int phaseIdx) const
Returns the average molar mass of the fluid phase.
Definition: porousmediumflow/2pnc/volumevariables.hh:348
SolidState solidState_
Definition: porousmediumflow/2pnc/volumevariables.hh:482
FluidState fluidState_
Definition: porousmediumflow/2pnc/volumevariables.hh:481
typename Traits::FluidSystem FluidSystem
Export fluid system type.
Definition: porousmediumflow/2pnc/volumevariables.hh:100
typename Traits::FluidState FluidState
Export fluid state type.
Definition: porousmediumflow/2pnc/volumevariables.hh:98
void completeFluidState(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv, FluidState &fluidState, SolidState &solidState)
Sets complete fluid state.
Definition: porousmediumflow/2pnc/volumevariables.hh:182
int wettingPhase() const
Returns the wetting phase index.
Definition: porousmediumflow/2pnc/volumevariables.hh:477
Scalar massFraction(int phaseIdx, int compIdx) const
Returns the mass fraction of a component in the phase.
Definition: porousmediumflow/2pnc/volumevariables.hh:462
Scalar pressure(int phaseIdx) const
Returns the effective pressure of a given phase within the control volume.
Definition: porousmediumflow/2pnc/volumevariables.hh:399
static constexpr bool useMoles()
Return whether moles or masses are balanced.
Definition: porousmediumflow/2pnc/volumevariables.hh:111
Scalar mobility(int phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition: porousmediumflow/2pnc/volumevariables.hh:418
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition: porousmediumflow/2pnc/volumevariables.hh:334
Definition: porousmediumflow/nonisothermal/volumevariables.hh:76
The isothermal base class.
Definition: porousmediumflow/volumevariables.hh:42
static constexpr int numFluidComponents()
Return number of components considered by the model.
Definition: porousmediumflow/volumevariables.hh:54
const PrimaryVariables & priVars() const
Returns the vector of primary variables.
Definition: porousmediumflow/volumevariables.hh:78
void update(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv)
Updates all quantities for a given control volume.
Definition: porousmediumflow/volumevariables.hh:66
Base class for the model specific class which provides access to all volume averaged quantities.
Base class for the model specific class which provides access to all volume averaged quantities.
The primary variable switch for the extended Richards model.