Commit 688168d0 authored by Alexander Kissinger's avatar Alexander Kissinger
Browse files

Commit of CO2 and CO2ni boxmodels, constitutive relationships and two

tests. (Reviewed by Lena)


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@9131 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 55040bb6
......@@ -82,6 +82,8 @@ AC_CONFIG_FILES([dumux.pc
test/boxmodels/2p2cni/Makefile
test/boxmodels/3p3c/Makefile
test/boxmodels/3p3cni/Makefile
test/boxmodels/co2/Makefile
test/boxmodels/co2ni/Makefile
test/boxmodels/mpnc/Makefile
test/boxmodels/richards/Makefile
test/common/Makefile
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2012 by *
* Institute for Modelling Hydraulic and Environmental Systems *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
#ifndef DUMUX_CO2_MODEL_HH
#define DUMUX_CO2_MODEL_HH
#include <dumux/boxmodels/2p2c/2p2cmodel.hh>
namespace Dumux
{
template<class TypeTag>
class CO2Model: public TwoPTwoCModel<TypeTag>
{
typedef TwoPTwoCModel<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, BaseModel) BaseType;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
enum {
numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
numComponents = GET_PROP_VALUE(TypeTag, NumComponents)
};
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum {
pressureIdx = Indices::pressureIdx,
switchIdx = Indices::switchIdx,
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
wCompIdx = Indices::wCompIdx,
nCompIdx = Indices::nCompIdx,
wPhaseOnly = Indices::wPhaseOnly,
nPhaseOnly = Indices::nPhaseOnly,
bothPhases = Indices::bothPhases,
pwSn = TwoPTwoCFormulation::pwSn,
pnSw = TwoPTwoCFormulation::pnSw,
formulation = GET_PROP_VALUE(TypeTag, Formulation)
};
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::ctype CoordScalar;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
public:
/*!
* \brief Update the static data of all vertices in the grid.
*
* \param curGlobalSol The current global solution
* \param oldGlobalSol The previous global solution
*/
void updateStaticData(SolutionVector &curGlobalSol,
const SolutionVector &oldGlobalSol)
{
bool wasSwitched = false;
for (unsigned i = 0; i < ParentType::staticVertexDat_.size(); ++i)
ParentType::staticVertexDat_[i].visited = false;
unsigned numDofs = this->numDofs();
unsigned numVertices = this->problem_().gridView().size(dim);
FVElementGeometry fvGeometry;
static VolumeVariables volVars;
ElementIterator eIt = this->gridView_().template begin<0> ();
const ElementIterator &eEndIt = this->gridView_().template end<0> ();
for (; eIt != eEndIt; ++eIt)
{
fvGeometry.update(this->gridView_(), *eIt);
for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
int globalIdx;
if (numDofs != numVertices)
globalIdx = this->elementMapper().map(*eIt);
else
globalIdx = this->vertexMapper().map(*eIt, scvIdx, dim);
if (ParentType::staticVertexDat_[globalIdx].visited)
continue;
ParentType::staticVertexDat_[globalIdx].visited = true;
volVars.update(curGlobalSol[globalIdx],
this->problem_(),
*eIt,
fvGeometry,
scvIdx,
false);
const GlobalPosition &globalPos = eIt->geometry().corner(scvIdx);
if (primaryVarSwitch_(curGlobalSol,
volVars,
globalIdx,
globalPos))
{
this->jacobianAssembler().markVertexRed(globalIdx);
wasSwitched = true;
}
}
}
// make sure that if there was a variable switch in an
// other partition we will also set the switch flag
// for our partition.
if (this->gridView_().comm().size() > 1)
wasSwitched = this->gridView_().comm().max(wasSwitched);
ParentType::setSwitched_(wasSwitched);
}
protected:
/*!
* \brief Set the old phase of all verts state to the current one.
*/
bool primaryVarSwitch_(SolutionVector &globalSol,
const VolumeVariables &volVars, int globalIdx,
const GlobalPosition &globalPos)
{
typename FluidSystem::ParameterCache paramCache;
// evaluate primary variable switch
bool wouldSwitch = false;
int phasePresence = ParentType::staticVertexDat_[globalIdx].phasePresence;
int newPhasePresence = phasePresence;
// check if a primary var switch is necessary
if (phasePresence == nPhaseOnly)
{
Scalar xnw = volVars.fluidState().moleFraction(nPhaseIdx, wCompIdx);
Scalar xnwMax = FluidSystem::equilibriumMoleFraction(volVars.fluidState(), paramCache, nPhaseIdx);
if(xnw > xnwMax)
wouldSwitch = true;
if (ParentType::staticVertexDat_[globalIdx].wasSwitched)
xnwMax *= 1.02;
//If mole fraction is higher than the equilibrium mole fraction make a phase switch
if(xnw > xnwMax)
{
// wetting phase appears
std::cout << "wetting phase appears at vertex " << globalIdx
<< ", coordinates: " << globalPos << ", xnw > xnwMax: "
<< xnw << " > "<< xnwMax << std::endl;
newPhasePresence = bothPhases;
if (formulation == pnSw)
globalSol[globalIdx][switchIdx] = 0.0;
else if (formulation == pwSn)
globalSol[globalIdx][switchIdx] = 1.0;
}
}
else if (phasePresence == wPhaseOnly)
{
Scalar xwn = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
Scalar xwnMax = FluidSystem::equilibriumMoleFraction(volVars.fluidState(), paramCache, wPhaseIdx);
//If mole fraction is higher than the equilibrium mole fraction make a phase switch
if(xwn > xwnMax)
wouldSwitch = true;
if (ParentType::staticVertexDat_[globalIdx].wasSwitched)
xwnMax *= 1.02;
if(xwn > xwnMax)
{
std::cout << "wetting phase appears at vertex " << globalIdx
<< ", coordinates: " << globalPos << ", xwn > xwnMax: "
<< xwn << " > "<< xwnMax << std::endl;
newPhasePresence = bothPhases;
if (formulation == pnSw)
globalSol[globalIdx][switchIdx] = 0.999;
else if (formulation == pwSn)
globalSol[globalIdx][switchIdx] = 0.001;
}
}
else if (phasePresence == bothPhases)
{
Scalar Smin = 0.0;
if (ParentType::staticVertexDat_[globalIdx].wasSwitched)
Smin = -0.01;
if (volVars.saturation(nPhaseIdx) <= Smin)
{
wouldSwitch = true;
// nonwetting phase disappears
std::cout << "Nonwetting phase disappears at vertex " << globalIdx
<< ", coordinates: " << globalPos << ", Sn: "
<< volVars.saturation(nPhaseIdx) << std::endl;
newPhasePresence = wPhaseOnly;
globalSol[globalIdx][switchIdx]
= volVars.fluidState().massFraction(wPhaseIdx, nCompIdx);
}
else if (volVars.saturation(wPhaseIdx) <= Smin)
{
wouldSwitch = true;
// wetting phase disappears
std::cout << "Wetting phase disappears at vertex " << globalIdx
<< ", coordinates: " << globalPos << ", Sw: "
<< volVars.saturation(wPhaseIdx) << std::endl;
newPhasePresence = nPhaseOnly;
globalSol[globalIdx][switchIdx]
= volVars.fluidState().massFraction(nPhaseIdx, wCompIdx);
}
}
ParentType::staticVertexDat_[globalIdx].phasePresence = newPhasePresence;
ParentType::staticVertexDat_[globalIdx].wasSwitched = wouldSwitch;
return phasePresence != newPhasePresence;
}
};
}
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2012 by *
* Institute for Modelling Hydraulic and Environmental Systems *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
#ifndef DUMUX_CO2_VOLUME_VARIABLES_HH
#define DUMUX_CO2_VOLUME_VARIABLES_HH
#include <dumux/boxmodels/2p2c/2p2cvolumevariables.hh>
namespace Dumux
{
template <class TypeTag>
class CO2VolumeVariables: public TwoPTwoCVolumeVariables<TypeTag>
{
typedef TwoPTwoCVolumeVariables<TypeTag> ParentType;
typedef BoxVolumeVariables<TypeTag> BaseClassType;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
// typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum {
numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
numComponents = GET_PROP_VALUE(TypeTag, NumComponents)
};
enum {
wCompIdx = Indices::wCompIdx,
nCompIdx = Indices::nCompIdx,
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx
};
// present phases
enum {
wPhaseOnly = Indices::lPhaseOnly,
nPhaseOnly = Indices::gPhaseOnly,
bothPhases = Indices::bothPhases
};
// formulations
enum {
formulation = GET_PROP_VALUE(TypeTag, Formulation),
pwSn = TwoPTwoCFormulation::pwSn,
pnSw = TwoPTwoCFormulation::pnSw
};
// primary variable indices
enum {
switchIdx = Indices::switchIdx,
pressureIdx = Indices::pressureIdx
};
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
enum { dim = GridView::dimension};
static const Scalar R; // universial nonwetting constant
public:
//! The type of the object returned by the fluidState() method
typedef Dumux::CompositionalFluidState<Scalar, FluidSystem> FluidState;
/*!
* \brief Update all quantities for a given control volume.
*
* \param priVars The primary variables
* \param problem The problem
* \param element The element
* \param fvGeometry The finite-volume geometry in the box scheme
* \param scvIdx The local index of the SCV (sub-control volume)
* \param isOldSol Evaluate function with solution of current or previous time step
*/
void update(const PrimaryVariables &priVars,
const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx,
const bool isOldSol)
{
// Update BoxVolVars but not 2p2cvolvars
// ToDo: Is BaseClassType the right name?
BaseClassType::update(priVars,
problem,
element,
fvGeometry,
scvIdx,
isOldSol);
unsigned numDofs = problem.model().numDofs();
unsigned numVertices = problem.gridView().size(dim);
int globalIdx;
if (numDofs != numVertices) // element data
globalIdx = problem.model().dofMapper().map(element);
else
globalIdx = problem.model().dofMapper().map(element, scvIdx, dim);
int phasePresence = problem.model().phasePresence(globalIdx, isOldSol);
Scalar temp = Implementation::temperature_(priVars, problem, element, fvGeometry, scvIdx);
ParentType::fluidState_.setTemperature(temp);
/////////////
// set the saturations
/////////////
Scalar Sn;
if (phasePresence == nPhaseOnly)
Sn = 1.0;
else if (phasePresence == wPhaseOnly) {
Sn = 0.0;
}
else if (phasePresence == bothPhases) {
if (formulation == pwSn)
Sn = priVars[switchIdx];
else if (formulation == pnSw)
Sn = 1.0 - priVars[switchIdx];
else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
}
else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
ParentType::fluidState_.setSaturation(wPhaseIdx, 1 - Sn);
ParentType::fluidState_.setSaturation(nPhaseIdx, Sn);
// capillary pressure parameters
const MaterialLawParams &materialParams =
problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
Scalar pC = MaterialLaw::pC(materialParams, 1 - Sn);
if (formulation == pwSn) {
ParentType::fluidState_.setPressure(wPhaseIdx, priVars[pressureIdx]);
ParentType::fluidState_.setPressure(nPhaseIdx, priVars[pressureIdx] + pC);
}
else if (formulation == pnSw) {
ParentType::fluidState_.setPressure(nPhaseIdx, priVars[pressureIdx]);
ParentType::fluidState_.setPressure(wPhaseIdx, priVars[pressureIdx] - pC);
}
else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
/////////////
// calculate the phase compositions
/////////////
typename FluidSystem::ParameterCache paramCache;
// calculate phase composition
if (phasePresence == bothPhases) {
//Get the equilibrium mole fractions from the FluidSystem and set them in the fluidState
//xCO2 = equilibrium mole fraction of CO2 in the liquid phase
//yH2O = equilibrium mole fraction of H2O in the gas phase
Scalar xwCO2 = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, wPhaseIdx);
Scalar xgH2O = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, nPhaseIdx);
Scalar xwH2O = 1 - xwCO2;
Scalar xgCO2 = 1 - xgH2O;
ParentType::fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xwH2O);
ParentType::fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwCO2);
ParentType::fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xgH2O);
ParentType::fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xgCO2);
//Get the phase densities from the FluidSystem and set them in the fluidState
Scalar rhoW = FluidSystem::density(ParentType::fluidState_, paramCache, wPhaseIdx);
Scalar rhoN = FluidSystem::density(ParentType::fluidState_, paramCache, nPhaseIdx);
ParentType::fluidState_.setDensity(wPhaseIdx, rhoW);
ParentType::fluidState_.setDensity(nPhaseIdx, rhoN);
//Get the phase viscosities from the FluidSystem and set them in the fluidState
Scalar muW = FluidSystem::viscosity(ParentType::fluidState_, paramCache, wPhaseIdx);
Scalar muN = FluidSystem::viscosity(ParentType::fluidState_, paramCache, nPhaseIdx);
ParentType::fluidState_.setViscosity(wPhaseIdx, muW);
ParentType::fluidState_.setViscosity(nPhaseIdx, muN);
}
else if (phasePresence == nPhaseOnly) {
// only the nonwetting phase is present, i.e. nonwetting phase
// composition is stored explicitly.
// extract _mass_ fractions in the nonwetting phase
Scalar massFractionN[numComponents];
massFractionN[wCompIdx] = priVars[switchIdx];
massFractionN[nCompIdx] = 1 - massFractionN[wCompIdx];
// calculate average molar mass of the nonwetting phase
Scalar M1 = FluidSystem::molarMass(wCompIdx);
Scalar M2 = FluidSystem::molarMass(nCompIdx);
Scalar X2 = massFractionN[nCompIdx];
Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
// convert mass to mole fractions and set the fluid state
ParentType::fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, massFractionN[wCompIdx]*avgMolarMass/M1);
ParentType::fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, massFractionN[nCompIdx]*avgMolarMass/M2);
// TODO give values for non-existing wetting phase
Scalar xwCO2 = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, wPhaseIdx);
Scalar xwH2O = 1 - xwCO2;
// Scalar xwCO2 = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, wPhaseIdx, nPhaseOnly);
// Scalar xwH2O = 1 - xwCO2;
ParentType::fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwCO2);
ParentType::fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xwH2O);
//Get the phase densities from the FluidSystem and set them in the fluidState
Scalar rhoW = FluidSystem::density(ParentType::fluidState_, paramCache, wPhaseIdx);
Scalar rhoN = FluidSystem::density(ParentType::fluidState_, paramCache, nPhaseIdx);
ParentType::fluidState_.setDensity(wPhaseIdx, rhoW);
ParentType::fluidState_.setDensity(nPhaseIdx, rhoN);
//Get the phase viscosities from the FluidSystem and set them in the fluidState
Scalar muW = FluidSystem::viscosity(ParentType::fluidState_, paramCache, wPhaseIdx);
Scalar muN = FluidSystem::viscosity(ParentType::fluidState_, paramCache, nPhaseIdx);
ParentType::fluidState_.setViscosity(wPhaseIdx, muW);
ParentType::fluidState_.setViscosity(nPhaseIdx, muN);
}
else if (phasePresence == wPhaseOnly) {
// only the wetting phase is present, i.e. wetting phase
// composition is stored explicitly.
// extract _mass_ fractions in the nonwetting phase
Scalar massFractionW[numComponents];
massFractionW[nCompIdx] = priVars[switchIdx];
massFractionW[wCompIdx] = 1 - massFractionW[nCompIdx];
// calculate average molar mass of the nonwetting phase
Scalar M1 = FluidSystem::molarMass(wCompIdx);
Scalar M2 = FluidSystem::molarMass(nCompIdx);
Scalar X2 = massFractionW[nCompIdx];
Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
// convert mass to mole fractions and set the fluid state
ParentType::fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, massFractionW[wCompIdx]*avgMolarMass/M1);
ParentType::fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, massFractionW[nCompIdx]*avgMolarMass/M2);
// TODO give values for non-existing nonwetting phase
Scalar xnH2O = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, nPhaseIdx);
Scalar xnCO2 = 1 - xnH2O; //FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, nPhaseIdx);
// Scalar xnH2O = FluidSystem::equilibriumMoleFraction(ParentType::fluidState_, paramCache, nPhaseIdx, wPhaseOnly);
// Scalar xnCO2 = 1 - xnH2O;
ParentType::fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnCO2);
ParentType::fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnH2O);
Scalar rhoW = FluidSystem::density(ParentType::fluidState_, paramCache, wPhaseIdx);
Scalar rhoN = FluidSystem::density(ParentType::fluidState_, paramCache, nPhaseIdx);
ParentType::fluidState_.setDensity(wPhaseIdx, rhoW);
ParentType::fluidState_.setDensity(nPhaseIdx, rhoN);
//Get the phase viscosities from the FluidSystem and set them in the fluidState
Scalar muW = FluidSystem::viscosity(ParentType::fluidState_, paramCache, wPhaseIdx);
Scalar muN = FluidSystem::viscosity(ParentType::fluidState_, paramCache, nPhaseIdx);
ParentType::fluidState_.setViscosity(wPhaseIdx, muW);
ParentType::fluidState_.setViscosity(nPhaseIdx, muN);
}
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
// compute and set the enthalpy
Scalar h = Implementation::enthalpy_(ParentType::fluidState_, paramCache, phaseIdx);
ParentType::fluidState_.setEnthalpy(phaseIdx, h);
}
paramCache.updateAll(ParentType::fluidState_);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
// relative permeabilities
Scalar kr;
if (phaseIdx == wPhaseIdx)
kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx));
else // ATTENTION: krn requires the liquid saturation
// as parameter!
kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx));
ParentType::relativePermeability_[phaseIdx] = kr;
Valgrind::CheckDefined(ParentType::relativePermeability_[phaseIdx]);
// binary diffusion coefficents
ParentType::diffCoeff_[phaseIdx] =
FluidSystem::binaryDiffusionCoefficient(ParentType::fluidState_,