diff --git a/dumux/material/constraintsolvers/compositionalflash.hh b/dumux/material/constraintsolvers/compositionalflash.hh
index 76fad086061f8ae9e3d7f3b30e363f798bdb8ce0..f5cb1e0af6cb526dba2c18f35309ef60c3ce4b3a 100644
--- a/dumux/material/constraintsolvers/compositionalflash.hh
+++ b/dumux/material/constraintsolvers/compositionalflash.hh
@@ -161,12 +161,12 @@ public:
* - Assign total concentration to the present phase
*
* \param fluidState The sequential fluid state
- * \param Z1 Feed mass fraction \f$\mathrm{[-]}\f$
+ * \param Z0 Feed mass fraction: Mass of first component per total mass \f$\mathrm{[-]}\f$
* \param phasePressure Vector holding the pressure \f$\mathrm{[Pa]}\f$
* \param presentPhaseIdx Subdomain Index = Indication which phase is present
* \param temperature Temperature \f$\mathrm{[K]}\f$
*/
- static void concentrationFlash1p2c(FluidState1p2c& fluidState, const Scalar& Z1,const Dune::FieldVector<Scalar,numPhases>
+ static void concentrationFlash1p2c(FluidState1p2c& fluidState, const Scalar& Z0,const Dune::FieldVector<Scalar,numPhases>
phasePressure,const int presentPhaseIdx, const Scalar& temperature)
{
// set the temperature, pressure
@@ -175,47 +175,15 @@ public:
fluidState.setPressure(phase1Idx, phasePressure[phase1Idx]);
fluidState.setPresentPhaseIdx(presentPhaseIdx);
- fluidState.setMassFraction(presentPhaseIdx,comp0Idx, Z1);
-
- // calculate mole fraction and average molar mass
- 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 == phase0Idx)
-// {
-//
-//// fluidState.setMassFraction(phase0Idx,comp0Idx, 0.;
-//
-//
-//
-//
-//
-//// fluidState.moleFractionWater_[phase1Idx] = 0.;
-//
-// fluidState.setPresentPhaseIdx(presentPhaseIdx);
-// }
-// else if (presentPhaseIdx == phase1Idx)
-// {
-// fluidState.setMassFraction[phase0Idx] = 0.;
-// fluidState.setMassFraction[phase1Idx] = Z1;
-//
-// // interested in nComp => 1-X1
-// 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_[phase1Idx] = 0.;
-//
-// fluidState.presentPhaseIdx_ = phase1Idx;
-// }
-// else
-// Dune::dgrave << __FILE__ <<": Twophase conditions in single-phase flash!"
-// << " Z1 is "<<Z1<< std::endl;
+ fluidState.setMassFraction(presentPhaseIdx,comp0Idx, Z0);
+
+ // transform mass to mole fractions
+ fluidState.setMoleFraction(presentPhaseIdx, comp0Idx, Z0 / FluidSystem::molarMass(comp0Idx)
+ / (Z0 / FluidSystem::molarMass(comp0Idx) + (1. - Z0) / FluidSystem::molarMass(comp1Idx)));
fluidState.setAverageMolarMass(presentPhaseIdx,
- fluidState.massFraction(presentPhaseIdx, comp0Idx)*FluidSystem::molarMass(comp0Idx)
- +fluidState.massFraction(presentPhaseIdx, comp1Idx)*FluidSystem::molarMass(comp1Idx));
+ fluidState.massFraction(presentPhaseIdx, comp0Idx) * FluidSystem::molarMass(comp0Idx)
+ + fluidState.massFraction(presentPhaseIdx, comp1Idx) * FluidSystem::molarMass(comp1Idx));
fluidState.setDensity(presentPhaseIdx, FluidSystem::density(fluidState, presentPhaseIdx));
fluidState.setMolarDensity(presentPhaseIdx, FluidSystem::molarDensity(fluidState, presentPhaseIdx));
@@ -251,17 +219,13 @@ public:
fluidState.setPressure(phase0Idx, phasePressure[phase0Idx]);
fluidState.setPressure(phase1Idx, phasePressure[phase1Idx]);
- //in contrast to the standard update() method, satflash() does not calculate nu.
- fluidState.setNu(phase0Idx, NAN);
- fluidState.setNu(phase1Idx, NAN);
-
//mole equilibrium ratios K for in case wPhase is reference phase
- double k1 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx); // = p^wComp_vap
- double k2 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx); // = H^nComp_w
+ double k00 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx);
+ double k01 = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx);
- // get mole fraction from equilibrium konstants
- fluidState.setMoleFraction(phase0Idx,comp0Idx, ((1. - k2) / (k1 -k2)));
- fluidState.setMoleFraction(phase1Idx,comp0Idx, (fluidState.moleFraction(phase0Idx,comp0Idx) * k1));
+ // get mole fraction from equilibrium constants
+ fluidState.setMoleFraction(phase0Idx,comp0Idx, ((1. - k01) / (k00 - k01)));
+ fluidState.setMoleFraction(phase1Idx,comp0Idx, (fluidState.moleFraction(phase0Idx,comp0Idx) * k00));
// transform mole to mass fractions
fluidState.setMassFraction(phase0Idx, comp0Idx,
diff --git a/dumux/porousmediumflow/2p2c/sequential/fvpressure.hh b/dumux/porousmediumflow/2p2c/sequential/fvpressure.hh
index 5f17f75288597aefd6528aeac1cc4e0f0a1b813c..177ad8714f2c432a669283590530878375c0bf78 100644
--- a/dumux/porousmediumflow/2p2c/sequential/fvpressure.hh
+++ b/dumux/porousmediumflow/2p2c/sequential/fvpressure.hh
@@ -919,23 +919,23 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
if(postTimeStep)
cellData.reset();
- // get the overall mass of component 1 Z1 = C^k / (C^1+C^2) [-]
- Scalar Z1 = cellData.massConcentration(wCompIdx)
+ // get the overall mass of first component Z0 = C^0 / (C^+C^1) [-]
+ Scalar Z0 = cellData.massConcentration(wCompIdx)
/ (cellData.massConcentration(wCompIdx)
+ cellData.massConcentration(nCompIdx));
// make shure only physical quantities enter flash calculation
- if(Z1 < 0. || Z1 > 1.)
+ if(Z0 < 0. || Z0 > 1.)
{
- Dune::dgrave << "Feed mass fraction unphysical: Z1 = " << Z1
+ Dune::dgrave << "Feed mass fraction unphysical: Z0 = " << Z0
<< " at global Idx " << eIdxGlobal
<< " , because totalConcentration(wCompIdx) = "
<< cellData.totalConcentration(wCompIdx)
<< " and totalConcentration(nCompIdx) = "
<< cellData.totalConcentration(nCompIdx)<< std::endl;
- if(Z1 < 0.)
+ if(Z0 < 0.)
{
- Z1 = 0.;
+ Z0 = 0.;
cellData.setTotalConcentration(wCompIdx, 0.);
problem().transportModel().totalConcentration(wCompIdx, eIdxGlobal) = 0.;
Dune::dgrave << "Regularize totalConcentration(wCompIdx) = "
@@ -943,7 +943,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
}
else
{
- Z1 = 1.;
+ Z0 = 1.;
cellData.setTotalConcentration(nCompIdx, 0.);
problem().transportModel().totalConcentration(nCompIdx,eIdxGlobal) = 0.;
Dune::dgrave << "Regularize totalConcentration(eIdxGlobal, nCompIdx) = "
@@ -973,7 +973,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
//complete fluid state
CompositionalFlash<Scalar, FluidSystem> flashSolver;
- flashSolver.concentrationFlash2p2c(fluidState, Z1, pressure, problem().spatialParams().porosity(elementI), temperature_);
+ flashSolver.concentrationFlash2p2c(fluidState, Z0, pressure, problem().spatialParams().porosity(elementI), temperature_);
// iterations part in case of enabled capillary pressure
Scalar pc(0.), oldPc(0.);
@@ -1006,7 +1006,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
//store old pc
oldPc = pc;
//update with better pressures
- flashSolver.concentrationFlash2p2c(fluidState, Z1, pressure,
+ flashSolver.concentrationFlash2p2c(fluidState, Z0, pressure,
problem().spatialParams().porosity(elementI), problem().temperatureAtPos(globalPos));
pc = MaterialLaw::pc(problem().spatialParams().materialLawParams(elementI),
fluidState.saturation(wPhaseIdx));
diff --git a/dumux/porousmediumflow/2p2c/sequential/fvpressurecompositional.hh b/dumux/porousmediumflow/2p2c/sequential/fvpressurecompositional.hh
index 0badb7ec78487a677ce6052172e3e98733b8bb5d..8f881ca7d034a19f7bd31e3c1646418343d0ada9 100644
--- a/dumux/porousmediumflow/2p2c/sequential/fvpressurecompositional.hh
+++ b/dumux/porousmediumflow/2p2c/sequential/fvpressurecompositional.hh
@@ -566,8 +566,8 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
}
else if (icFormulation == Indices::concentration) // concentration initial condition
{
- Scalar Z1_0 = problem_.initConcentration(element);
- flashSolver.concentrationFlash2p2c(fluidState, Z1_0, pressure,
+ Scalar Z0 = problem_.initConcentration(element);
+ flashSolver.concentrationFlash2p2c(fluidState, Z0, pressure,
problem_.spatialParams().porosity(element), temperature_);
}
}
@@ -607,7 +607,7 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
}
else if (icFormulation == Indices::concentration) // concentration initial condition
{
- Scalar Z1_0 = problem_.initConcentration(element);
+ Scalar Z0 = problem_.initConcentration(element);
// If total concentrations are given at the boundary, saturation is unknown.
// This may affect pc and hence p_alpha and hence again saturation -> iteration.
@@ -641,7 +641,7 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
//store old pc
Scalar oldPc = pc;
//update with better pressures
- flashSolver.concentrationFlash2p2c(fluidState, Z1_0, pressure,
+ flashSolver.concentrationFlash2p2c(fluidState, Z0, pressure,
problem_.spatialParams().porosity(element), problem_.temperatureAtPos(globalPos));
pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element),
fluidState.saturation(wPhaseIdx));
@@ -659,7 +659,7 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
{
pressure[wPhaseIdx] = pressure[nPhaseIdx]
= this->pressure()[eIdxGlobal];
- flashSolver.concentrationFlash2p2c(fluidState, Z1_0,
+ flashSolver.concentrationFlash2p2c(fluidState, Z0,
pressure, problem_.spatialParams().porosity(element), temperature_);
}
} //end conc initial condition
@@ -803,8 +803,8 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
// numerical derivative of fluid volume with respect to pressure
PhaseVector p_(incp);
p_ += pressure;
- Scalar Z1 = mass[0] / mass.one_norm();
- flashSolver.concentrationFlash2p2c(updFluidState, Z1,
+ Scalar Z0 = mass[0] / mass.one_norm();
+ flashSolver.concentrationFlash2p2c(updFluidState, Z0,
p_, problem_.spatialParams().porosity(element), temperature_);
specificVolume=0.; // = \sum_{\alpha} \nu_{\alpha} / \rho_{\alpha}
@@ -818,7 +818,7 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
Dune::dinfo << "dv_dp larger 0 at Idx " << eIdxGlobal << " , try and invert secant"<< std::endl;
p_ -= 2*incp;
- flashSolver.concentrationFlash2p2c(updFluidState, Z1,
+ flashSolver.concentrationFlash2p2c(updFluidState, Z0,
p_, problem_.spatialParams().porosity(element), temperature_);
specificVolume=0.; // = \sum_{\alpha} \nu_{\alpha} / \rho_{\alpha}
@@ -831,7 +831,7 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
{
Dune::dwarn << "dv_dp still larger 0 after inverting secant at idx"<< eIdxGlobal<< std::endl;
p_ += 2*incp;
- flashSolver.concentrationFlash2p2c(updFluidState, Z1,
+ flashSolver.concentrationFlash2p2c(updFluidState, Z0,
p_, problem_.spatialParams().porosity(element), temperature_);
// neglect effects of phase split, only regard changes in phase densities
specificVolume=0.; // = \sum_{\alpha} \nu_{\alpha} / \rho_{\alpha}
@@ -851,9 +851,9 @@ void FVPressureCompositional<TypeTag>::volumeDerivatives(const GlobalPosition& g
for (int compIdx = 0; compIdx<numComponents; compIdx++)
{
mass[compIdx] += massIncrement[compIdx];
- Z1 = mass[0] / mass.one_norm();
+ Z0 = mass[0] / mass.one_norm();
- flashSolver.concentrationFlash2p2c(updFluidState, Z1,
+ flashSolver.concentrationFlash2p2c(updFluidState, Z0,
pressure, problem_.spatialParams().porosity(element), temperature_);
specificVolume=0.; // = \sum_{\alpha} \nu_{\alpha} / \rho_{\alpha}
diff --git a/dumux/porousmediumflow/2p2c/sequential/fvpressuremultiphysics.hh b/dumux/porousmediumflow/2p2c/sequential/fvpressuremultiphysics.hh
index 81c92a3acff159dd91246c5d46a7b5954edcd0f5..83ec75b1e1688d97920a6617e88bc5f0832e81bd 100644
--- a/dumux/porousmediumflow/2p2c/sequential/fvpressuremultiphysics.hh
+++ b/dumux/porousmediumflow/2p2c/sequential/fvpressuremultiphysics.hh
@@ -425,11 +425,11 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pStorage(Dune::FieldVector<Scalar,
p_[wPhaseIdx] += cellDataI.pressure(wPhaseIdx);
Scalar sumC = (cellDataI.massConcentration(wCompIdx) + cellDataI.massConcentration(nCompIdx));
- Scalar Z1 = cellDataI.massConcentration(wCompIdx) / sumC;
+ Scalar Z0 = cellDataI.massConcentration(wCompIdx) / sumC;
// initialize simple fluidstate object
PseudoOnePTwoCFluidState<Scalar, FluidSystem> pseudoFluidState;
CompositionalFlash<Scalar, FluidSystem> flashSolver;
- flashSolver.concentrationFlash1p2c(pseudoFluidState, Z1, p_, cellDataI.subdomain(),
+ flashSolver.concentrationFlash1p2c(pseudoFluidState, Z0, p_, cellDataI.subdomain(),
cellDataI.temperature(wPhaseIdx));
Scalar v_ = 1. / pseudoFluidState.density(presentPhaseIdx);
@@ -441,7 +441,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pStorage(Dune::FieldVector<Scalar,
Dune::dinfo << "dv_dp larger 0 at Idx " << eIdxGlobalI << " , try and invert secant"<< std::endl;
p_ -= 2*incp;
- flashSolver.concentrationFlash1p2c(pseudoFluidState, Z1, p_, cellDataI.subdomain(),
+ flashSolver.concentrationFlash1p2c(pseudoFluidState, Z0, p_, cellDataI.subdomain(),
cellDataI.temperature(wPhaseIdx));
v_ = 1. / pseudoFluidState.density(presentPhaseIdx);
cellDataI.dv_dp() = (sumC * ( v_ - (1. /cellDataI.density(presentPhaseIdx)))) /incp;
@@ -950,13 +950,13 @@ void FVPressure2P2CMultiPhysics<TypeTag>::update1pMaterialLawsInElement(const El
pressure[nPhaseIdx] = this->pressure(eIdxGlobal);
}
- // get the overall mass of component 1: Z1 = C^k / (C^1+C^2) [-]
+ // get the overall mass of first component: Z0 = C^0 / (C^0+C^1) [-]
Scalar sumConc = cellData.massConcentration(wCompIdx)
+ cellData.massConcentration(nCompIdx);
- Scalar Z1 = cellData.massConcentration(wCompIdx)/ sumConc;
+ Scalar Z0 = cellData.massConcentration(wCompIdx)/ sumConc;
CompositionalFlash<Scalar, FluidSystem> flashSolver;
- flashSolver.concentrationFlash1p2c(pseudoFluidState, Z1, pressure, presentPhaseIdx, problem().temperatureAtPos(globalPos));
+ flashSolver.concentrationFlash1p2c(pseudoFluidState, Z0, pressure, presentPhaseIdx, problem().temperatureAtPos(globalPos));
// write stuff in fluidstate
assert(presentPhaseIdx == pseudoFluidState.presentPhaseIdx());
diff --git a/dumux/porousmediumflow/2p2c/sequential/fvtransport.hh b/dumux/porousmediumflow/2p2c/sequential/fvtransport.hh
index 8e713b11f1ece326e3938795924aff7e47ba9dc1..0ad6215e1a6002c4672592ff665140b834612e9f 100644
--- a/dumux/porousmediumflow/2p2c/sequential/fvtransport.hh
+++ b/dumux/porousmediumflow/2p2c/sequential/fvtransport.hh
@@ -1153,8 +1153,8 @@ void FVTransport2P2C<TypeTag>::evalBoundary(GlobalPosition globalPosFace,
{
// saturation and hence pc and hence corresponding pressure unknown
pressBound[wPhaseIdx] = pressBound[nPhaseIdx] = primaryVariablesOnBoundary[Indices::pressureEqIdx];
- Scalar Z1Bound = primaryVariablesOnBoundary[contiWEqIdx];
- flashSolver.concentrationFlash2p2c(BCfluidState, Z1Bound, pressBound,
+ Scalar Z0Bound = primaryVariablesOnBoundary[contiWEqIdx];
+ flashSolver.concentrationFlash2p2c(BCfluidState, Z0Bound, pressBound,
problem().spatialParams().porosity(element), problem().temperatureAtPos(globalPosFace));
if(GET_PROP_VALUE(TypeTag, EnableCapillarity))
@@ -1187,7 +1187,7 @@ void FVTransport2P2C<TypeTag>::evalBoundary(GlobalPosition globalPosFace,
//store old pc
Scalar oldPc = pcBound;
//update with better pressures
- flashSolver.concentrationFlash2p2c(BCfluidState, Z1Bound, pressBound,
+ flashSolver.concentrationFlash2p2c(BCfluidState, Z0Bound, pressBound,
problem().spatialParams().porosity(element), problem().temperatureAtPos(globalPosFace));
pcBound = MaterialLaw::pc(problem().spatialParams().materialLawParams(element),
BCfluidState.saturation(wPhaseIdx));