Commit 08e525da authored by Timo Koch's avatar Timo Koch
Browse files

[2p][1pnc][richardsnc] Improve naming / definition of phases and components and their indices

Fluid systems exported phase indices for wetting and nonwetting components. The fluid system
is not the right place to decide about wettability. This introduces wettingPhase() in the spatial
params which returns the index of the phase in the fluidsystem we consider the wetting phase.
This way wettability can change temporally and spatially.
Indices throughout the 2p/2p2c/2pnc nc and mpnc where renamed accordingly.
parent 9dca5357
......@@ -128,7 +128,7 @@ NEW_PROP_TAG(EvaluatePermeabilityAtScvfIP);
// Additional properties used by the 2pnc and 2pncmin models:
//////////////////////////////////////////////////////////////
NEW_PROP_TAG(Chemistry); //!< The chemistry class with which solves equlibrium reactions
NEW_PROP_TAG(SetMoleFractionsForWettingPhase); //!< Set the mole fraction in the wetting or non-wetting phase
NEW_PROP_TAG(SetMoleFractionsForFirstPhase); //!< Set the mole fraction in the wetting or non-wetting phase
//////////////////////////////////////////////////////////////
// Additional properties used by the richards model
......
......@@ -41,9 +41,12 @@ template<class Scalar, class ThermalConductivityModel, class FluidSystem>
class PlotThermalConductivityModel
{
using FluidState = CompositionalFluidState<Scalar, FluidSystem>;
enum {
wPhaseIdx = FluidSystem::wPhaseIdx,
nPhaseIdx = FluidSystem::nPhaseIdx
// phase indices
enum
{
phase0Idx = FluidSystem::phase0Idx,
phase1Idx = FluidSystem::phase1Idx
};
public:
......@@ -61,10 +64,10 @@ public:
{
FluidState fluidstate;
fluidstate.setTemperature(temperature);
fluidstate.setPressure(wPhaseIdx, pressure);
fluidstate.setPressure(nPhaseIdx, pressure);
lambdaW_ = FluidSystem::thermalConductivity(fluidstate, wPhaseIdx);
lambdaN_ = FluidSystem::thermalConductivity(fluidstate, nPhaseIdx);
fluidstate.setPressure(phase0Idx, pressure);
fluidstate.setPressure(phase1Idx, pressure);
lambdaW_ = FluidSystem::thermalConductivity(fluidstate, phase0Idx);
lambdaN_ = FluidSystem::thermalConductivity(fluidstate, phase1Idx);
}
/*!
......
......@@ -26,7 +26,6 @@
#include <cmath>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/exceptions.hh>
#include <dumux/discretization/methods.hh>
......@@ -46,37 +45,36 @@ enum ElectroChemistryModel { Ochs, Acosta };
* \brief This class calculates source terms and current densities for fuel cells
* with the electrochemical models suggested by Ochs (2008) \cite ochs2008 or Acosta et al. (2006) \cite A3:acosta:2006
* \todo TODO: Scalar type should be extracted from VolumeVariables!
* \todo TODO: This shouldn't depend on grid and discretization!!
*/
template <class Scalar, class Indices, class FluidSystem, class FVGridGeometry, ElectroChemistryModel electroChemistryModel>
class ElectroChemistry
{
using FVElementGeometry = typename FVGridGeometry::LocalView;
using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using Constant = Dumux::Constants<Scalar>;
enum {
//indices of the components
wCompIdx = FluidSystem::wCompIdx, //major component of the liquid phase
nCompIdx = FluidSystem::nCompIdx, //major component of the gas phase
O2Idx = wCompIdx + 2
};
enum {
enum
{
//indices of the primary variables
pressureIdx = Indices::pressureIdx, //gas-phase pressure
switchIdx = Indices::switchIdx, //liquid saturation or mole fraction
pressureIdx = Indices::pressureIdx, // gas-phase pressure
switchIdx = Indices::switchIdx, // liquid saturation or mole fraction
};
enum {
enum
{
//equation indices
conti0EqIdx = Indices::conti0EqIdx,
contiH2OEqIdx = conti0EqIdx + wCompIdx,
contiO2EqIdx = conti0EqIdx + wCompIdx + 2,
contiH2OEqIdx = Indices::conti0EqIdx + FluidSystem::H2OIdx,
contiO2EqIdx = Indices::conti0EqIdx + FluidSystem::O2Idx,
};
static constexpr bool isBox = FVGridGeometry::discMethod == DiscretizationMethod::box;
enum { dofCodim = isBox ? GridView::dimension : 0 };
enum
{
// phase indices
liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
gasPhaseIdx = FluidSystem::gasPhaseIdx
};
using GridView = typename FVGridGeometry::GridView;
static constexpr bool isBox = FVGridGeometry::discMethod == DiscretizationMethod::box;
using GlobalPosition = typename Dune::FieldVector<typename GridView::ctype, GridView::dimensionworld>;
public:
......@@ -92,9 +90,9 @@ public:
static void reactionSource(SourceValues &values, Scalar currentDensity,
const std::string& paramGroup = "")
{
//correction to account for actually relevant reaction area
//current density has to be devided by the half length of the box
//\todo Do we have multiply with the electrochemically active surface area (ECSA) here instead?
// correction to account for actually relevant reaction area
// current density has to be devided by the half length of the box
// \todo Do we have multiply with the electrochemically active surface area (ECSA) here instead?
static Scalar gridYMax = getParamFromGroup<GlobalPosition>(paramGroup, "Grid.UpperRight")[1];
static Scalar nCellsY = getParamFromGroup<GlobalPosition>(paramGroup, "Grid.Cells")[1];
......@@ -143,7 +141,7 @@ public:
Scalar activationLossesNext = calculateActivationLosses_(volVars, currentDensity+deltaCurrentDensity);
Scalar concentrationLosses = calculateConcentrationLosses_(volVars);
Scalar activationLossesDiff = activationLossesNext - activationLosses;
Scalar sw = volVars.saturation(FluidSystem::wPhaseIdx);
Scalar sw = volVars.saturation(liquidPhaseIdx);
if(electroChemistryModel == ElectroChemistryModel::Acosta)
{
......@@ -198,11 +196,11 @@ private:
static Scalar transferCoefficient = getParam<Scalar>("ElectroChemistry.TransferCoefficient");
//Saturation sw for Acosta calculation
Scalar sw = volVars.saturation(FluidSystem::wPhaseIdx);
Scalar sw = volVars.saturation(liquidPhaseIdx);
//Calculate prefactor
Scalar preFactor = Constant::R*volVars.fluidState().temperature()/transferCoefficient/Constant::F/numElectrons;
//Get partial pressure of O2 in the gas phase
Scalar pO2 = volVars.pressure(FluidSystem::nPhaseIdx) * volVars.fluidState().moleFraction(FluidSystem::nPhaseIdx, O2Idx);
Scalar pO2 = volVars.pressure(gasPhaseIdx) * volVars.fluidState().moleFraction(gasPhaseIdx, FluidSystem::O2Idx);
Scalar losses = 0.0;
//Calculate activation losses
......@@ -241,7 +239,7 @@ private:
//Calculate preFactor
Scalar preFactor = Constant::R*volVars.temperature()/transferCoefficient/Constant::F/numElectrons;
//Get partial pressure of O2 in the gas phase
Scalar pO2 = volVars.pressure(FluidSystem::nPhaseIdx) * volVars.fluidState().moleFraction(FluidSystem::nPhaseIdx, O2Idx);
Scalar pO2 = volVars.pressure(gasPhaseIdx) * volVars.fluidState().moleFraction(gasPhaseIdx, FluidSystem::O2Idx);
Scalar losses = 0.0;
//Calculate concentration losses
......
......@@ -24,7 +24,6 @@
#ifndef DUMUX_ELECTROCHEMISTRY_NI_HH
#define DUMUX_ELECTROCHEMISTRY_NI_HH
#include <dumux/common/properties.hh>
#include <dumux/material/constants.hh>
#include <dumux/material/chemistry/electrochemistry/electrochemistry.hh>
......@@ -36,30 +35,23 @@ namespace Dumux {
* with the electrochemical models suggested by Ochs (2008) \cite ochs2008 or Acosta (2006) \cite A3:acosta:2006
* for the non-isothermal case.
* \todo TODO: Scalar type should be extracted from VolumeVariables!
* \todo TODO: This shouldn't depend on discretization and grid!!
*/
template <class Scalar, class Indices, class FVGridGeometry, ElectroChemistryModel electroChemistryModel>
class ElectroChemistryNI : public ElectroChemistry<Scalar, Indices, FVGridGeometry, electroChemistryModel>
{
using ParentType = ElectroChemistry<Scalar, Indices, FVGridGeometry, electroChemistryModel>;
using GridView = typename FVGridGeometry::GridView;
using Constant = Constants<Scalar>;
enum {
//indices of the components
wCompIdx = Indices::wCompIdx, //major component of the liquid phase
nCompIdx = Indices::nCompIdx, //major component of the gas phase
O2Idx = wCompIdx + 2
};
enum { //equation indices
conti0EqIdx = Indices::conti0EqIdx,
contiH2OEqIdx = conti0EqIdx + wCompIdx,
contiO2EqIdx = conti0EqIdx + wCompIdx + 2,
energyEqIdx = Indices::energyEqIdx, //energy equation
//equation indices
contiH2OEqIdx = Indices::conti0EqIdx + FluidSystem::H2OIdx,
contiO2EqIdx = Indices::conti0EqIdx + FluidSystem::O2Idx,
energyEqIdx = Indices::energyEqIdx, //energy equation
};
using GridView = typename FVGridGeometry::GridView;
static constexpr bool isBox = FVGridGeometry::discMethod == DiscretizationMethod::box;
enum { dofCodim = isBox ? GridView::dimension : 0 };
using GlobalPosition = typename Dune::FieldVector<typename GridView::ctype, GridView::dimensionworld>;
using CellVector = typename Dune::FieldVector<typename GridView::ctype, GridView::dimension>;
......
......@@ -47,10 +47,10 @@ class CompositionalFlash
};
enum{
wPhaseIdx = FluidSystem::wPhaseIdx,
nPhaseIdx = FluidSystem::nPhaseIdx,
wCompIdx = FluidSystem::wCompIdx,
nCompIdx = FluidSystem::nCompIdx
phase0Idx = FluidSystem::phase0Idx,
phase1Idx = FluidSystem::phase1Idx,
comp0Idx = FluidSystem::comp0Idx,
comp1Idx = FluidSystem::comp1Idx
};
public:
......@@ -82,93 +82,93 @@ public:
{
// set the temperature, pressure
fluidState.setTemperature(temperature);
fluidState.setPressure(wPhaseIdx, phasePressure[wPhaseIdx]);
fluidState.setPressure(nPhaseIdx, phasePressure[nPhaseIdx]);
fluidState.setPressure(phase0Idx, phasePressure[phase0Idx]);
fluidState.setPressure(phase1Idx, phasePressure[phase1Idx]);
//mole equilibrium ratios K for in case wPhase is reference phase
double k1 = FluidSystem::fugacityCoefficient(fluidState, wPhaseIdx, wCompIdx); // = p^wComp_vap
double k2 = FluidSystem::fugacityCoefficient(fluidState, wPhaseIdx, nCompIdx); // = H^nComp_w
double k1 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx); // = p^wComp_vap
double k2 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx); // = H^nComp_w
// get mole fraction from equilibrium konstants
fluidState.setMoleFraction(wPhaseIdx,wCompIdx, ((1. - k2) / (k1 -k2)));
fluidState.setMoleFraction(nPhaseIdx,wCompIdx, (fluidState.moleFraction(wPhaseIdx,wCompIdx) * k1));
fluidState.setMoleFraction(phase0Idx,comp0Idx, ((1. - k2) / (k1 -k2)));
fluidState.setMoleFraction(phase1Idx,comp0Idx, (fluidState.moleFraction(phase0Idx,comp0Idx) * k1));
// transform mole to mass fractions
fluidState.setMassFraction(wPhaseIdx, wCompIdx,
(fluidState.moleFraction(wPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
/ ( fluidState.moleFraction(wPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
+ (1.-fluidState.moleFraction(wPhaseIdx,wCompIdx)) * FluidSystem::molarMass(nCompIdx) )));
fluidState.setMassFraction(nPhaseIdx,wCompIdx,
(fluidState.moleFraction(nPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
/ ( fluidState.moleFraction(nPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
+ (1.-fluidState.moleFraction(nPhaseIdx,wCompIdx)) * FluidSystem::molarMass(nCompIdx) )));
fluidState.setMassFraction(phase0Idx, comp0Idx,
(fluidState.moleFraction(phase0Idx,comp0Idx) * FluidSystem::molarMass(comp0Idx)
/ ( fluidState.moleFraction(phase0Idx,comp0Idx) * FluidSystem::molarMass(comp0Idx)
+ (1.-fluidState.moleFraction(phase0Idx,comp0Idx)) * FluidSystem::molarMass(comp1Idx) )));
fluidState.setMassFraction(phase1Idx,comp0Idx,
(fluidState.moleFraction(phase1Idx,comp0Idx) * FluidSystem::molarMass(comp0Idx)
/ ( fluidState.moleFraction(phase1Idx,comp0Idx) * FluidSystem::molarMass(comp0Idx)
+ (1.-fluidState.moleFraction(phase1Idx,comp0Idx)) * FluidSystem::molarMass(comp1Idx) )));
//mass equilibrium ratios
Scalar equilRatio_[numPhases][numComponents];
equilRatio_[nPhaseIdx][wCompIdx] = fluidState.massFraction(nPhaseIdx,wCompIdx)
/ fluidState.massFraction(wPhaseIdx, wCompIdx); // = Xn1 / Xw1 = K1
equilRatio_[nPhaseIdx][nCompIdx] = (1.-fluidState.massFraction(nPhaseIdx, wCompIdx))
/ (1.-fluidState.massFraction(wPhaseIdx, wCompIdx)); // =(1.-Xn1) / (1.-Xw1) = K2
equilRatio_[wPhaseIdx][nCompIdx] = equilRatio_[wPhaseIdx][wCompIdx] = 1.;
equilRatio_[phase1Idx][comp0Idx] = fluidState.massFraction(phase1Idx,comp0Idx)
/ fluidState.massFraction(phase0Idx, comp0Idx); // = Xn1 / Xw1 = K1
equilRatio_[phase1Idx][comp1Idx] = (1.-fluidState.massFraction(phase1Idx, comp0Idx))
/ (1.-fluidState.massFraction(phase0Idx, comp0Idx)); // =(1.-Xn1) / (1.-Xw1) = K2
equilRatio_[phase0Idx][comp1Idx] = equilRatio_[phase0Idx][comp0Idx] = 1.;
// phase fraction of nPhase [mass/totalmass]
fluidState.setNu(nPhaseIdx, 0.);
fluidState.setNu(phase1Idx, 0.);
// check if there is enough of component 1 to form a phase
if (Z1 > fluidState.massFraction(nPhaseIdx,wCompIdx)
&& Z1 < fluidState.massFraction(wPhaseIdx,wCompIdx))
fluidState.setNu(nPhaseIdx, -((equilRatio_[nPhaseIdx][wCompIdx]-1)*Z1 + (equilRatio_[nPhaseIdx][nCompIdx]-1)*(1-Z1))
/ (equilRatio_[nPhaseIdx][wCompIdx]-1) / (equilRatio_[nPhaseIdx][nCompIdx] -1));
else if (Z1 <= fluidState.massFraction(nPhaseIdx,wCompIdx)) // too little wComp to form a phase
if (Z1 > fluidState.massFraction(phase1Idx,comp0Idx)
&& Z1 < fluidState.massFraction(phase0Idx,comp0Idx))
fluidState.setNu(phase1Idx, -((equilRatio_[phase1Idx][comp0Idx]-1)*Z1 + (equilRatio_[phase1Idx][comp1Idx]-1)*(1-Z1))
/ (equilRatio_[phase1Idx][comp0Idx]-1) / (equilRatio_[phase1Idx][comp1Idx] -1));
else if (Z1 <= fluidState.massFraction(phase1Idx,comp0Idx)) // too little wComp to form a phase
{
fluidState.setNu(nPhaseIdx, 1.); // only nPhase
fluidState.setMassFraction(nPhaseIdx,wCompIdx, Z1); // hence, assign complete mass dissolved into nPhase
fluidState.setNu(phase1Idx, 1.); // only nPhase
fluidState.setMassFraction(phase1Idx,comp0Idx, Z1); // hence, assign complete mass dissolved into nPhase
// store as moleFractions
Scalar xw_n = Z1 /*=Xw_n*/ / FluidSystem::molarMass(wCompIdx); // = moles of compIdx
xw_n /= ( Z1 /*=Xw_n*/ / FluidSystem::molarMass(wCompIdx)
+(1- Z1 /*=Xn_n*/) / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
Scalar xw_n = Z1 /*=Xw_n*/ / FluidSystem::molarMass(comp0Idx); // = moles of compIdx
xw_n /= ( Z1 /*=Xw_n*/ / FluidSystem::molarMass(comp0Idx)
+(1- Z1 /*=Xn_n*/) / FluidSystem::molarMass(comp1Idx) ); // /= total moles in phase
fluidState.setMoleFraction(nPhaseIdx,wCompIdx, xw_n);
fluidState.setMoleFraction(phase1Idx,comp0Idx, xw_n);
// // opposing non-present phase is already set to equilibrium mass fraction
// fluidState.setMassFraction(wPhaseIdx,wCompIdx, 1.); // non present phase is set to be pure
// fluidState.setMoleFraction(wPhaseIdx,wCompIdx, 1.); // non present phase is set to be pure
// fluidState.setMassFraction(phase0Idx,comp0Idx, 1.); // non present phase is set to be pure
// fluidState.setMoleFraction(phase0Idx,comp0Idx, 1.); // non present phase is set to be pure
}
else // (Z1 >= Xw1) => no nPhase
{
fluidState.setNu(nPhaseIdx, 0.); // no second phase
fluidState.setMassFraction(wPhaseIdx, wCompIdx, Z1);
fluidState.setNu(phase1Idx, 0.); // no second phase
fluidState.setMassFraction(phase0Idx, comp0Idx, Z1);
// store as moleFractions
Scalar xw_w = Z1 /*=Xw_w*/ / FluidSystem::molarMass(wCompIdx); // = moles of compIdx
xw_w /= ( Z1 /*=Xw_w*/ / FluidSystem::molarMass(wCompIdx)
+(1- Z1 /*=Xn_w*/) / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
fluidState.setMoleFraction(wPhaseIdx, wCompIdx, xw_w);
Scalar xw_w = Z1 /*=Xw_w*/ / FluidSystem::molarMass(comp0Idx); // = moles of compIdx
xw_w /= ( Z1 /*=Xw_w*/ / FluidSystem::molarMass(comp0Idx)
+(1- Z1 /*=Xn_w*/) / FluidSystem::molarMass(comp1Idx) ); // /= total moles in phase
fluidState.setMoleFraction(phase0Idx, comp0Idx, xw_w);
// // opposing non-present phase is already set to equilibrium mass fraction
// fluidState.setMassFraction(nPhaseIdx,wCompIdx, 0.); // non present phase is set to be pure
// fluidState.setMoleFraction(nPhaseIdx,wCompIdx, 0.); // non present phase is set to be pure
// fluidState.setMassFraction(phase1Idx,comp0Idx, 0.); // non present phase is set to be pure
// fluidState.setMoleFraction(phase1Idx,comp0Idx, 0.); // non present phase is set to be pure
}
// complete array of mass fractions
fluidState.setMassFraction(wPhaseIdx, nCompIdx, 1. - fluidState.massFraction(wPhaseIdx,wCompIdx));
fluidState.setMassFraction(nPhaseIdx, nCompIdx, 1. - fluidState.massFraction(nPhaseIdx,wCompIdx));
fluidState.setMassFraction(phase0Idx, comp1Idx, 1. - fluidState.massFraction(phase0Idx,comp0Idx));
fluidState.setMassFraction(phase1Idx, comp1Idx, 1. - fluidState.massFraction(phase1Idx,comp0Idx));
// complete array of mole fractions
fluidState.setMoleFraction(wPhaseIdx, nCompIdx, 1. - fluidState.moleFraction(wPhaseIdx,wCompIdx));
fluidState.setMoleFraction(nPhaseIdx, nCompIdx, 1. - fluidState.moleFraction(nPhaseIdx,wCompIdx));
fluidState.setMoleFraction(phase0Idx, comp1Idx, 1. - fluidState.moleFraction(phase0Idx,comp0Idx));
fluidState.setMoleFraction(phase1Idx, comp1Idx, 1. - fluidState.moleFraction(phase1Idx,comp0Idx));
// complete phase mass fractions
fluidState.setNu(wPhaseIdx, 1. - fluidState.phaseMassFraction(nPhaseIdx));
fluidState.setNu(phase0Idx, 1. - fluidState.phaseMassFraction(phase1Idx));
// get densities with correct composition
fluidState.setDensity(wPhaseIdx, FluidSystem::density(fluidState, wPhaseIdx));
fluidState.setDensity(nPhaseIdx, FluidSystem::density(fluidState, nPhaseIdx));
fluidState.setDensity(phase0Idx, FluidSystem::density(fluidState, phase0Idx));
fluidState.setDensity(phase1Idx, FluidSystem::density(fluidState, phase1Idx));
Scalar sw = fluidState.phaseMassFraction(wPhaseIdx) / fluidState.density(wPhaseIdx);
sw /= (fluidState.phaseMassFraction(wPhaseIdx)/fluidState.density(wPhaseIdx)
+ fluidState.phaseMassFraction(nPhaseIdx)/fluidState.density(nPhaseIdx));
fluidState.setSaturation(wPhaseIdx, sw);
Scalar sw = fluidState.phaseMassFraction(phase0Idx) / fluidState.density(phase0Idx);
sw /= (fluidState.phaseMassFraction(phase0Idx)/fluidState.density(phase0Idx)
+ fluidState.phaseMassFraction(phase1Idx)/fluidState.density(phase1Idx));
fluidState.setSaturation(phase0Idx, sw);
}
/*! The simplest possible update routine for 1p2c "flash" calculations
......@@ -188,51 +188,51 @@ public:
{
// set the temperature, pressure
fluidState.setTemperature(temperature);
fluidState.setPressure(wPhaseIdx, phasePressure[wPhaseIdx]);
fluidState.setPressure(nPhaseIdx, phasePressure[nPhaseIdx]);
fluidState.setPressure(phase0Idx, phasePressure[phase0Idx]);
fluidState.setPressure(phase1Idx, phasePressure[phase1Idx]);
fluidState.setPresentPhaseIdx(presentPhaseIdx);
fluidState.setMassFraction(presentPhaseIdx,wCompIdx, Z1);
fluidState.setMassFraction(presentPhaseIdx,comp0Idx, Z1);
// calculate mole fraction and average molar mass
Scalar xw_alpha= Z1 / FluidSystem::molarMass(wCompIdx);
xw_alpha /= ( Z1 / FluidSystem::molarMass(wCompIdx)
+ (1.-Z1) / FluidSystem::molarMass(nCompIdx));
fluidState.setMoleFraction(presentPhaseIdx, wCompIdx, xw_alpha);
Scalar xw_alpha= Z1 / FluidSystem::molarMass(comp0Idx);
xw_alpha /= ( Z1 / FluidSystem::molarMass(comp0Idx)
+ (1.-Z1) / FluidSystem::molarMass(comp1Idx));
fluidState.setMoleFraction(presentPhaseIdx, comp0Idx, xw_alpha);
// if (presentPhaseIdx == wPhaseIdx)
// if (presentPhaseIdx == phase0Idx)
// {
//
//// fluidState.setMassFraction(wPhaseIdx,wCompIdx, 0.;
//// fluidState.setMassFraction(phase0Idx,comp0Idx, 0.;
//
//
//
//
//
//// fluidState.moleFractionWater_[nPhaseIdx] = 0.;
//// fluidState.moleFractionWater_[phase1Idx] = 0.;
//
// fluidState.setPresentPhaseIdx(presentPhaseIdx);
// }
// else if (presentPhaseIdx == nPhaseIdx)
// else if (presentPhaseIdx == phase1Idx)
// {
// fluidState.setMassFraction[wPhaseIdx] = 0.;
// fluidState.setMassFraction[nPhaseIdx] = Z1;
// fluidState.setMassFraction[phase0Idx] = 0.;
// fluidState.setMassFraction[phase1Idx] = Z1;
//
// // interested in nComp => 1-X1
// fluidState.moleFractionWater_[nPhaseIdx] = ( Z1 / FluidSystem::molarMass(0) ); // = moles of compIdx
// fluidState.moleFractionWater_[nPhaseIdx] /= (Z1/ FluidSystem::molarMass(0)
// fluidState.moleFractionWater_[phase1Idx] = ( Z1 / FluidSystem::molarMass(0) ); // = moles of compIdx
// fluidState.moleFractionWater_[phase1Idx] /= (Z1/ FluidSystem::molarMass(0)
// + (1.-Z1) / FluidSystem::molarMass(1) ); // /= total moles in phase
// fluidState.moleFractionWater_[nPhaseIdx] = 0.;
// fluidState.moleFractionWater_[phase1Idx] = 0.;
//
// fluidState.presentPhaseIdx_ = nPhaseIdx;
// fluidState.presentPhaseIdx_ = phase1Idx;
// }
// else
// Dune::dgrave << __FILE__ <<": Twophase conditions in single-phase flash!"
// << " Z1 is "<<Z1<< std::endl;
fluidState.setAverageMolarMass(presentPhaseIdx,
fluidState.massFraction(presentPhaseIdx, wCompIdx)*FluidSystem::molarMass(wCompIdx)
+fluidState.massFraction(presentPhaseIdx, nCompIdx)*FluidSystem::molarMass(nCompIdx));
fluidState.massFraction(presentPhaseIdx, comp0Idx)*FluidSystem::molarMass(comp0Idx)
+fluidState.massFraction(presentPhaseIdx, comp1Idx)*FluidSystem::molarMass(comp1Idx));
fluidState.setDensity(presentPhaseIdx, FluidSystem::density(fluidState, presentPhaseIdx));
}
......@@ -263,44 +263,44 @@ public:
{
// set the temperature, pressure
fluidState.setTemperature(temperature);
fluidState.setPressure(wPhaseIdx, phasePressure[wPhaseIdx]);
fluidState.setPressure(nPhaseIdx, phasePressure[nPhaseIdx]);
fluidState.setPressure(phase0Idx, phasePressure[phase0Idx]);
fluidState.setPressure(phase1Idx, phasePressure[phase1Idx]);
//in contrast to the standard update() method, satflash() does not calculate nu.
fluidState.setNu(wPhaseIdx, NAN);
fluidState.setNu(nPhaseIdx, NAN);
fluidState.setNu(phase0Idx, NAN);
fluidState.setNu(phase1Idx, NAN);
//mole equilibrium ratios K for in case wPhase is reference phase
double k1 = FluidSystem::fugacityCoefficient(fluidState, wPhaseIdx, wCompIdx); // = p^wComp_vap
double k2 = FluidSystem::fugacityCoefficient(fluidState, wPhaseIdx, nCompIdx); // = H^nComp_w
double k1 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx); // = p^wComp_vap
double k2 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx); // = H^nComp_w
// get mole fraction from equilibrium konstants
fluidState.setMoleFraction(wPhaseIdx,wCompIdx, ((1. - k2) / (k1 -k2)));
fluidState.setMoleFraction(nPhaseIdx,wCompIdx, (fluidState.moleFraction(wPhaseIdx,wCompIdx) * k1));
fluidState.setMoleFraction(phase0Idx,comp0Idx, ((1. - k2) / (k1 -k2)));
fluidState.setMoleFraction(phase1Idx,comp0Idx, (fluidState.moleFraction(phase0Idx,comp0Idx) * k1));
// transform mole to mass fractions
fluidState.setMassFraction(wPhaseIdx, wCompIdx,
(fluidState.moleFraction(wPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
/ ( fluidState.moleFraction(wPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
+ (1.-fluidState.moleFraction(wPhaseIdx,wCompIdx)) * FluidSystem::molarMass(nCompIdx) )));
fluidState.setMassFraction(nPhaseIdx,wCompIdx,
(fluidState.moleFraction(nPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
/ ( fluidState.moleFraction(nPhaseIdx,wCompIdx) * FluidSystem::molarMass(wCompIdx)
+ (1.-fluidState.moleFraction(nPhaseIdx,wCompIdx)) * FluidSystem::molarMass(nCompIdx) )));
fluidState.setMassFraction(phase0Idx, comp0Idx,
(fluidState.moleFraction(phase0Idx,comp0Idx) * FluidSystem::molarMass(comp0Idx)
/ ( fluidState.moleFraction(phase0Idx,comp0Idx) * FluidSystem::molarMass(comp0Idx)
+ (1.-fluidState.moleFraction(phase0Idx,comp0Idx)) * FluidSystem::molarMass(comp1Idx) )));
fluidState.setMassFraction(phase1Idx,comp0Idx,
(fluidState.moleFraction(phase1Idx,comp0Idx) * FluidSystem::molarMass(comp0Idx)
/ ( fluidState.moleFraction(phase1Idx,comp0Idx) * FluidSystem::molarMass(comp0Idx)
+ (1.-fluidState.moleFraction(phase1Idx,comp0Idx)) * FluidSystem::molarMass(comp1Idx) )));
// complete array of mass fractions
fluidState.setMassFraction(wPhaseIdx, nCompIdx, 1. - fluidState.massFraction(wPhaseIdx,wCompIdx));
fluidState.setMassFraction(nPhaseIdx, nCompIdx, 1. - fluidState.massFraction(nPhaseIdx,wCompIdx));
fluidState.setMassFraction(phase0Idx, comp1Idx, 1. - fluidState.massFraction(phase0Idx,comp0Idx));
fluidState.setMassFraction(phase1Idx, comp1Idx, 1. - fluidState.massFraction(phase1Idx,comp0Idx));
// complete array of mole fractions
fluidState.setMoleFraction(wPhaseIdx, nCompIdx, 1. - fluidState.moleFraction(wPhaseIdx,wCompIdx));
fluidState.setMoleFraction(nPhaseIdx, nCompIdx, 1. - fluidState.moleFraction(nPhaseIdx,wCompIdx));
fluidState.setMoleFraction(phase0Idx, comp1Idx, 1. - fluidState.moleFraction(phase0Idx,comp0Idx));
fluidState.setMoleFraction(phase1Idx, comp1Idx, 1. - fluidState.moleFraction(phase1Idx,comp0Idx));
// get densities with correct composition
fluidState.setDensity(wPhaseIdx, FluidSystem::density(fluidState, wPhaseIdx));
fluidState.setDensity(nPhaseIdx, FluidSystem::density(fluidState, nPhaseIdx));
fluidState.setDensity(phase0Idx, FluidSystem::density(fluidState, phase0Idx));
fluidState.setDensity(phase1Idx, FluidSystem::density(fluidState, phase1Idx));
// set saturation
fluidState.setSaturation(wPhaseIdx, saturation);
fluidState.setSaturation(phase0Idx, saturation);
}
//@}
};
......
......@@ -57,7 +57,7 @@ namespace Dumux {
* - if the setViscosity parameter is true, also dynamic viscosities of *all* phases
* - if the setEnthalpy parameter is true, also specific enthalpies of *all* phases
*/
template <class Scalar, class FluidSystem, bool useKelvinEquation = false>
template <class Scalar, class FluidSystem>
class MiscibleMultiPhaseComposition
{
static constexpr int numPhases = FluidSystem::numPhases;
......@@ -140,9 +140,9 @@ public:
}
}
//set the additional equations for the numComponents-numMajorComponents
// set the additional equations for the numComponents-numMajorComponents
// this is only relevant if numphases != numcomponents e.g. in a 2pnc system
//Components, of which the molefractions are known, set to molefraction(knownCompIdx)=xKnown
// Components, of which the molefractions are known, set to molefraction(knownCompIdx)=xKnown
for(int knownCompIdx = 0; knownCompIdx < numComponents-numMajorComponents; ++knownCompIdx)
{
int rowIdx = numComponents + numPhases + knownCompIdx;
......@@ -156,70 +156,31 @@ public:
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
int col1Idx = compIdx;
Scalar entryPhase0 = 0.0;
for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
Scalar entry = fluidState.fugacityCoefficient(phaseIdx, compIdx)
* fluidState.pressure(phaseIdx);
// modify the saturation vapor pressure of the wetting component by the Kelvin equation
if (compIdx == FluidSystem::wCompIdx
&& phaseIdx == FluidSystem::wPhaseIdx
&& useKelvinEquation)
{
// a new fluidState is needed, because mole fractions are unknown
FluidState purePhaseFluidState;
// assign all phase pressures, needed for capillary pressure
for (unsigned int idx = 0; idx < numPhases; ++idx)
{
purePhaseFluidState.setPressure(idx, fluidState.pressure(idx));
}
purePhaseFluidState.setTemperature(fluidState.temperature());
purePhaseFluidState.setMoleFraction(phaseIdx, compIdx, 1.0);
entry = FluidSystem::kelvinVaporPressure(purePhaseFluidState, phaseIdx, compIdx);
}
const auto entryPhase0 = fluidState.fugacityCoefficient(0, compIdx)*fluidState.pressure(0);
if (phaseIdx == 0)
{
entryPhase0 = entry;
}
else
{
int rowIdx = (phaseIdx - 1)*numComponents + compIdx;
int col2Idx = phaseIdx*numComponents + compIdx;