Commit 6732681f authored by Benjamin Faigle's avatar Benjamin Faigle
Browse files

corrected fluidstate to changes in new fluidsystems solubility method

redone multiphyiscs module to work with new fluidsystem

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7158 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 8ced7dfc
......@@ -27,7 +27,6 @@
#ifndef DUMUX_DEC2P2C_FLUID_STATE_HH
#define DUMUX_DEC2P2C_FLUID_STATE_HH
#include <dumux/material/old_fluidsystems/fluidstate.hh>
#include <dumux/decoupled/2p2c/2p2cproperties.hh>
namespace Dumux
......@@ -42,8 +41,7 @@ namespace Dumux
* \tparam TypeTag The property Type Tag
*/
template <class TypeTag>
class DecoupledTwoPTwoCFluidState : public FluidState<typename GET_PROP_TYPE(TypeTag,
PTAG(Scalar)), DecoupledTwoPTwoCFluidState<TypeTag> >
class DecoupledTwoPTwoCFluidState
{
typedef DecoupledTwoPTwoCFluidState<TypeTag> ThisType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
......@@ -97,50 +95,64 @@ public:
//mole equilibrium ratios K for in case wPhase is reference phase
double k1 = FluidSystem::fugacityCoefficient(*this, wPhaseIdx, wCompIdx)
/ phasePressure_[nPhaseIdx]; // = p^wComp_vap / p_n
double k2 = FluidSystem::fugacityCoefficient(*this, wPhaseIdx, nCompIdx)
/ phasePressure_[nPhaseIdx]; // = H^nComp_w / p_n
double k1 = FluidSystem::fugacityCoefficient(*this, wPhaseIdx, wCompIdx); // = p^wComp_vap
double k2 = FluidSystem::fugacityCoefficient(*this, wPhaseIdx, nCompIdx); // = H^nComp_w
// get mole fraction from equilibrium konstants
Scalar xw1 = (1. - k2) / (k1 -k2);
Scalar xn1 = xw1 * k1;
moleFraction_[wPhaseIdx][wCompIdx] = (1. - k2) / (k1 -k2);
moleFraction_[nPhaseIdx][wCompIdx] = moleFraction_[wPhaseIdx][wCompIdx] * k1;
// transform mole to mass fractions
massfrac_[wPhaseIdx][wCompIdx] = xw1 * FluidSystem::molarMass(wCompIdx)
/ ( xw1 * FluidSystem::molarMass(wCompIdx) + (1.-xw1) * FluidSystem::molarMass(nCompIdx) );
massfrac_[nPhaseIdx][wCompIdx] = xn1 * FluidSystem::molarMass(wCompIdx)
/ ( xn1 * FluidSystem::molarMass(wCompIdx) + (1.-xn1) * FluidSystem::molarMass(nCompIdx) );
massFraction_[wPhaseIdx][wCompIdx] = moleFraction_[wPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx)
/ ( moleFraction_[wPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx) + (1.-moleFraction_[wPhaseIdx][wCompIdx]) * FluidSystem::molarMass(nCompIdx) );
massFraction_[nPhaseIdx][wCompIdx] = moleFraction_[nPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx)
/ ( moleFraction_[nPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx) + (1.-moleFraction_[nPhaseIdx][wCompIdx]) * FluidSystem::molarMass(nCompIdx) );
//mass equilibrium ratios
equilRatio_[nPhaseIdx][wCompIdx] = massfrac_[nPhaseIdx][wCompIdx] / massfrac_[wPhaseIdx][wCompIdx]; // = Xn1 / Xw1 = K1
equilRatio_[nPhaseIdx][nCompIdx] = (1.-massfrac_[nPhaseIdx][wCompIdx])/ (1.-massfrac_[wPhaseIdx][wCompIdx]); // =(1.-Xn1) / (1.-Xw1) = K2
equilRatio_[nPhaseIdx][wCompIdx] = massFraction_[nPhaseIdx][wCompIdx] / massFraction_[wPhaseIdx][wCompIdx]; // = Xn1 / Xw1 = K1
equilRatio_[nPhaseIdx][nCompIdx] = (1.-massFraction_[nPhaseIdx][wCompIdx])/ (1.-massFraction_[wPhaseIdx][wCompIdx]); // =(1.-Xn1) / (1.-Xw1) = K2
equilRatio_[wPhaseIdx][nCompIdx] = equilRatio_[wPhaseIdx][wCompIdx] = 1.;
// phase fraction of nPhase [mass/totalmass]
nu_[nPhaseIdx] = 0;
// check if there is enough of component 1 to form a phase
if (Z1 > massfrac_[nPhaseIdx][wCompIdx] && Z1 < massfrac_[wPhaseIdx][wCompIdx])
if (Z1 > massFraction_[nPhaseIdx][wCompIdx] && Z1 < massFraction_[wPhaseIdx][wCompIdx])
nu_[nPhaseIdx] = -((equilRatio_[nPhaseIdx][wCompIdx]-1)*Z1 + (equilRatio_[nPhaseIdx][nCompIdx]-1)*(1-Z1)) / (equilRatio_[nPhaseIdx][wCompIdx]-1) / (equilRatio_[nPhaseIdx][nCompIdx] -1);
else if (Z1 <= massfrac_[nPhaseIdx][wCompIdx]) // too little wComp to form a phase
else if (Z1 <= massFraction_[nPhaseIdx][wCompIdx]) // too little wComp to form a phase
{
nu_[nPhaseIdx] = 1; // only nPhase
massfrac_[nPhaseIdx][wCompIdx] = Z1; // hence, assign complete mass soluted into nPhase
massfrac_[wPhaseIdx][wCompIdx] = 1.;
massFraction_[nPhaseIdx][wCompIdx] = Z1; // hence, assign complete mass soluted into nPhase
// store as moleFractions
moleFraction_[nPhaseIdx][wCompIdx] = ( massFraction_[nPhaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx) ); // = moles of compIdx
moleFraction_[nPhaseIdx][wCompIdx] /= ( massFraction_[nPhaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx)
+ massFraction_[nPhaseIdx][nCompIdx] / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
// corresponding phase
massFraction_[wPhaseIdx][wCompIdx] = 1.;
moleFraction_[wPhaseIdx][wCompIdx] = 1.;
}
else // (Z1 >= Xw1) => no nPhase
{
nu_[nPhaseIdx] = 0; // no second phase
massfrac_[wPhaseIdx][wCompIdx] = Z1;
massfrac_[nPhaseIdx][wCompIdx] = 0.;
massFraction_[wPhaseIdx][wCompIdx] = Z1;
// store as moleFractions
moleFraction_[wPhaseIdx][wCompIdx] = ( massFraction_[wPhaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx) ); // = moles of compIdx
moleFraction_[wPhaseIdx][wCompIdx] /= ( massFraction_[wPhaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx)
+ massFraction_[wPhaseIdx][nCompIdx] / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
massFraction_[nPhaseIdx][wCompIdx] = 0.;
moleFraction_[nPhaseIdx][wCompIdx] = 0.;
}
// complete array of mass fractions
massfrac_[wPhaseIdx][nCompIdx] = 1. - massfrac_[wPhaseIdx][wCompIdx];
massfrac_[nPhaseIdx][nCompIdx] = 1. - massfrac_[nPhaseIdx][wCompIdx];
massFraction_[wPhaseIdx][nCompIdx] = 1. - massFraction_[wPhaseIdx][wCompIdx];
massFraction_[nPhaseIdx][nCompIdx] = 1. - massFraction_[nPhaseIdx][wCompIdx];
// complete array of mole fractions
moleFraction_[wPhaseIdx][nCompIdx] = 1. - moleFraction_[wPhaseIdx][wCompIdx];
moleFraction_[nPhaseIdx][nCompIdx] = 1. - moleFraction_[nPhaseIdx][wCompIdx];
// complete phase mass fractions
nu_[wPhaseIdx] = 1. - nu_[nPhaseIdx];
......@@ -153,11 +165,11 @@ public:
Sw_ /= (nu_[wPhaseIdx]/density_[wPhaseIdx] + nu_[nPhaseIdx]/density_[nPhaseIdx]);
massConcentration_[wCompIdx] =
poro * (massfrac_[wPhaseIdx][wCompIdx] * Sw_ * density_[wPhaseIdx]
+ massfrac_[nPhaseIdx][wCompIdx] * (1.-Sw_) * density_[nPhaseIdx]);
poro * (massFraction_[wPhaseIdx][wCompIdx] * Sw_ * density_[wPhaseIdx]
+ massFraction_[nPhaseIdx][wCompIdx] * (1.-Sw_) * density_[nPhaseIdx]);
massConcentration_[nCompIdx] =
poro * (massfrac_[wPhaseIdx][nCompIdx] * Sw_ * density_[wPhaseIdx]
+ massfrac_[nPhaseIdx][nCompIdx] * (1-Sw_) * density_[nPhaseIdx]);
poro * (massFraction_[wPhaseIdx][nCompIdx] * Sw_ * density_[wPhaseIdx]
+ massFraction_[nPhaseIdx][nCompIdx] * (1-Sw_) * density_[nPhaseIdx]);
}
//! a flash routine for 2p2c if the saturation instead of total concentration is known.
......@@ -196,23 +208,26 @@ public:
/ phasePressure_[nPhaseIdx];
// get mole fraction from equilibrium konstants
Scalar xw1 = (1. - k2) / (k1 -k2);
Scalar xn1 = xw1 * k1;
moleFraction_[wPhaseIdx][wCompIdx] = (1. - k2) / (k1 -k2);
moleFraction_[nPhaseIdx][wCompIdx] = moleFraction_[wPhaseIdx][wCompIdx] * k1;
// transform mole to mass fractions
massfrac_[wPhaseIdx][wCompIdx] = xw1 * FluidSystem::molarMass(wCompIdx)
/ ( xw1 * FluidSystem::molarMass(wCompIdx) + (1.-xw1) * FluidSystem::molarMass(nCompIdx) );
massfrac_[nPhaseIdx][wCompIdx] = xn1 * FluidSystem::molarMass(wCompIdx)
/ ( xn1 * FluidSystem::molarMass(wCompIdx) + (1.-xn1) * FluidSystem::molarMass(nCompIdx) );
massFraction_[wPhaseIdx][wCompIdx] = moleFraction_[wPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx)
/ ( moleFraction_[wPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx) + (1.-moleFraction_[wPhaseIdx][wCompIdx]) * FluidSystem::molarMass(nCompIdx) );
massFraction_[nPhaseIdx][wCompIdx] = moleFraction_[nPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx)
/ ( moleFraction_[nPhaseIdx][wCompIdx] * FluidSystem::molarMass(wCompIdx) + (1.-moleFraction_[nPhaseIdx][wCompIdx]) * FluidSystem::molarMass(nCompIdx) );
// complete array of mass fractions
massfrac_[wPhaseIdx][nCompIdx] = 1. - massfrac_[wPhaseIdx][wCompIdx];
massfrac_[nPhaseIdx][nCompIdx] = 1. - massfrac_[nPhaseIdx][wCompIdx];
massFraction_[wPhaseIdx][nCompIdx] = 1. - massFraction_[wPhaseIdx][wCompIdx];
massFraction_[nPhaseIdx][nCompIdx] = 1. - massFraction_[nPhaseIdx][wCompIdx];
// complete array of mole fractions
moleFraction_[wPhaseIdx][nCompIdx] = 1. - moleFraction_[wPhaseIdx][wCompIdx];
moleFraction_[nPhaseIdx][nCompIdx] = 1. - moleFraction_[nPhaseIdx][wCompIdx];
//mass equilibrium ratios
equilRatio_[nPhaseIdx][wCompIdx] = massfrac_[nPhaseIdx][wCompIdx] / massfrac_[wPhaseIdx][wCompIdx]; // = Xn1 / Xw1 = K1
equilRatio_[nPhaseIdx][nCompIdx] = (1.-massfrac_[nPhaseIdx][wCompIdx])/ (1.-massfrac_[wPhaseIdx][wCompIdx]); // =(1.-Xn1) / (1.-Xw1) = K2
equilRatio_[nPhaseIdx][wCompIdx] = massFraction_[nPhaseIdx][wCompIdx] / massFraction_[wPhaseIdx][wCompIdx]; // = Xn1 / Xw1 = K1
equilRatio_[nPhaseIdx][nCompIdx] = (1.-massFraction_[nPhaseIdx][wCompIdx])/ (1.-massFraction_[wPhaseIdx][wCompIdx]); // =(1.-Xn1) / (1.-Xw1) = K2
equilRatio_[wPhaseIdx][nCompIdx] = equilRatio_[wPhaseIdx][wCompIdx] = 1.;
// get densities with correct composition
......@@ -220,11 +235,11 @@ public:
density_[nPhaseIdx] = FluidSystem::density(*this, nPhaseIdx);
massConcentration_[wCompIdx] =
poro * (massfrac_[wPhaseIdx][wCompIdx] * Sw_ * density_[wPhaseIdx]
+ massfrac_[nPhaseIdx][wCompIdx] * (1.-Sw_) * density_[nPhaseIdx]);
poro * (massFraction_[wPhaseIdx][wCompIdx] * Sw_ * density_[wPhaseIdx]
+ massFraction_[nPhaseIdx][wCompIdx] * (1.-Sw_) * density_[nPhaseIdx]);
massConcentration_[nCompIdx] =
poro * (massfrac_[wPhaseIdx][nCompIdx] * Sw_ * density_[wPhaseIdx]
+ massfrac_[nPhaseIdx][nCompIdx] * (1-Sw_) * density_[nPhaseIdx]);
poro * (massFraction_[wPhaseIdx][nCompIdx] * Sw_ * density_[wPhaseIdx]
+ massFraction_[nPhaseIdx][nCompIdx] * (1-Sw_) * density_[nPhaseIdx]);
}
//@}
/*!
......@@ -252,7 +267,7 @@ public:
*/
Scalar massFraction(int phaseIdx, int compIdx) const
{
return massfrac_[phaseIdx][compIdx];
return massFraction_[phaseIdx][compIdx];
}
......@@ -264,11 +279,7 @@ public:
*/
Scalar moleFraction(int phaseIdx, int compIdx) const
{
// as the moass fractions are calculated, it is used to determine the mole fractions
double moleFrac_ = ( massfrac_[phaseIdx][compIdx] / FluidSystem::molarMass(compIdx) ); // = moles of compIdx
moleFrac_ /= ( massfrac_[phaseIdx][wCompIdx] / FluidSystem::molarMass(wCompIdx)
+ massfrac_[phaseIdx][nCompIdx] / FluidSystem::molarMass(nCompIdx) ); // /= total moles in phase
return moleFrac_;
return moleFraction_[phaseIdx][compIdx];
}
/*!
......@@ -334,6 +345,24 @@ public:
Scalar temperature(int phaseIdx) const
{ return temperature_; };
/*!
* \brief Return the average molar mass of a phase.
*
* This is the sum of all molar masses times their respective mole
* fractions in the phase.
*
* Unit: \f$\mathrm{[kg/m^3]}\f$
*/
Scalar averageMolarMass(int phaseIdx) const
{
Scalar averageMolarMass = 0;
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
averageMolarMass += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx);
}
return averageMolarMass;
}
/*!
* \brief Returns the phase mass fraction. phase mass per total mass \f$\mathrm{[kg/kg]}\f$.
*
......@@ -353,14 +382,17 @@ public:
//@}
private:
Scalar density_[numPhases];
Scalar massConcentration_[numComponents];
Scalar massfrac_[numPhases][numComponents];
Scalar equilRatio_[numPhases][numComponents];
Scalar phasePressure_[numPhases];
Scalar temperature_;
Scalar Sw_;
Scalar nu_[numPhases]; //phase mass fraction
Scalar density_[numPhases];
Scalar massFraction_[numPhases][numComponents];
Scalar moleFraction_[numPhases][numComponents];
Scalar equilRatio_[numPhases][numComponents];
Scalar averageMolarMass_[numPhases];
};
} // end namepace
......
......@@ -51,7 +51,7 @@ namespace Dumux
* Isothermal conditions and local thermodynamic
* equilibrium are assumed. Gravity is included.
* \f[
c_{total}\frac{\partial p}{\partial t} + \sum_{\kappa} \frac{\partial v_{total}}{\partial C^{\kappa}} \nabla \cdot \left( \sum_{\alpha} X^{\kappa}_{\alpha} \varrho_{alpha} \bf{v}_{\alpha}\right)
c_{total}\frac{\partial p}{\partial t} + \sum_{\kappa} \frac{\partial v_{total}}{\partial C^{\kappa}} \nabla \cdot \left( \sum_{\alpha} X^{\kappa}_{\alpha} \varrho_{\alpha} \bf{v}_{\alpha}\right)
= \sum_{\kappa} \frac{\partial v_{total}}{\partial C^{\kappa}} q^{\kappa},
* \f]
* where \f$\bf{v}_{\alpha} = - \lambda_{\alpha} \bf{K} \left(\nabla p_{\alpha} + \rho_{\alpha} \bf{g} \right) \f$.
......@@ -786,7 +786,6 @@ void FVPressure2P2C<TypeTag>::assemble(bool first)
// create a fluid state for the boundary
FluidState BCfluidState;
Scalar temperatureBC = problem().temperatureAtPos(globalPosFace);
//read boundary values
PrimaryVariables primaryVariablesOnBoundary(NAN);
problem().dirichlet(primaryVariablesOnBoundary, *isIt);
......
......@@ -606,7 +606,6 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
//read boundary values
PrimaryVariables primaryVariablesOnBoundary(NAN);
problem().dirichlet(primaryVariablesOnBoundary, *isIt);
Scalar temperatureBC = problem().temperatureAtPos(globalPosFace);
if (first)
{
......@@ -634,13 +633,13 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
Scalar lambdaNWBound = 0.;
Scalar densityWBound =
FluidSystem::phaseDensity(wPhaseIdx, temperatureBC, pressBC[wPhaseIdx], BCfluidState);
Scalar densityNWBound =
FluidSystem::phaseDensity(nPhaseIdx, temperatureBC, pressBC[nPhaseIdx], BCfluidState);
Scalar viscosityWBound =
FluidSystem::phaseViscosity(wPhaseIdx, temperatureBC, pressBC[wPhaseIdx], BCfluidState);
Scalar viscosityNWBound =
FluidSystem::phaseViscosity(nPhaseIdx, temperatureBC, pressBC[nPhaseIdx], BCfluidState);
FluidSystem::density(BCfluidState, wPhaseIdx);
Scalar densityNWBound =
FluidSystem::density(BCfluidState, nPhaseIdx);
Scalar viscosityWBound =
FluidSystem::viscosity(BCfluidState, wPhaseIdx);
Scalar viscosityNWBound =
FluidSystem::viscosity(BCfluidState, nPhaseIdx);
// mobility at the boundary
switch (GET_PROP_VALUE(TypeTag, PTAG(BoundaryMobility)))
......@@ -756,7 +755,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
else
{
densityW = densityWBound;
dV_w = (dv_dC1 * BCfluidState.massFrac(wPhaseIdx, wCompIdx) + dv_dC2 * BCfluidState.massFrac(wPhaseIdx, nCompIdx));
dV_w = (dv_dC1 * BCfluidState.massFraction(wPhaseIdx, wCompIdx) + dv_dC2 * BCfluidState.massFraction(wPhaseIdx, nCompIdx));
dV_w *= densityW;
lambdaW = lambdaWBound;
}
......@@ -771,7 +770,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
else
{
densityNW = densityNWBound;
dV_n = (dv_dC1 * BCfluidState.massFrac(nPhaseIdx, wCompIdx) + dv_dC2 * BCfluidState.massFrac(nPhaseIdx, nCompIdx));
dV_n = (dv_dC1 * BCfluidState.massFraction(nPhaseIdx, wCompIdx) + dv_dC2 * BCfluidState.massFraction(nPhaseIdx, nCompIdx));
dV_n *= densityNW;
lambdaNW = lambdaNWBound;
}
......@@ -857,12 +856,10 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
/ (problem().variables().totalConcentration(globalIdxI, wCompIdx)
+ problem().variables().totalConcentration(globalIdxI, nCompIdx));
PseudoOnePTwoCFluidState<TypeTag> pseudoFluidState;
pseudoFluidState.update(Z1, p_, problem().variables().saturation(globalIdxI));
pseudoFluidState.update(Z1, p_, problem().variables().saturation(globalIdxI), problem().temperatureAtPos(globalPos));
if(problem().variables().saturation(globalIdxI)==1) //only w-phase
{
Scalar v_w_ = 1. / FluidSystem::phaseDensity(wPhaseIdx,
problem().temperatureAtPos(globalPos),
p_, pseudoFluidState);
Scalar v_w_ = 1. / FluidSystem::density(pseudoFluidState, wPhaseIdx);
problem().variables().dv_dp(globalIdxI) = (problem().variables().totalConcentration(globalIdxI, wCompIdx)
+ problem().variables().totalConcentration(globalIdxI, nCompIdx))
* ( v_w_ - 1./densityWI) / incp;
......@@ -870,9 +867,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
if (problem().variables().dv_dp(globalIdxI) > 0) // this is not physically possible!!
{
p_ -= 2*incp;
Scalar v_w_ = 1. / FluidSystem::phaseDensity(wPhaseIdx,
problem().temperatureAtPos(globalPos),
p_, pseudoFluidState);
Scalar v_w_ = 1. / FluidSystem::density(pseudoFluidState, wPhaseIdx);
problem().variables().dv_dp(globalIdxI) = (problem().variables().totalConcentration(globalIdxI, wCompIdx)
+ problem().variables().totalConcentration(globalIdxI, nCompIdx))
* ( v_w_ - 1./densityWI) / -incp;
......@@ -880,9 +875,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
}
if(problem().variables().saturation(globalIdxI)==0.) //only nw-phase
{
Scalar v_n_ = 1. / FluidSystem::phaseDensity(nPhaseIdx,
problem().temperatureAtPos(globalPos),
p_, pseudoFluidState);
Scalar v_n_ = 1. / FluidSystem::density(pseudoFluidState, nPhaseIdx);
problem().variables().dv_dp(globalIdxI) = (problem().variables().totalConcentration(globalIdxI, wCompIdx)
+ problem().variables().totalConcentration(globalIdxI, nCompIdx))
* ( v_n_ - 1./densityNWI) / incp;
......@@ -890,9 +883,7 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
if (problem().variables().dv_dp(globalIdxI) > 0) // this is not physically possible!!
{
p_ -= 2*incp;
Scalar v_n_ = 1. / FluidSystem::phaseDensity(nPhaseIdx,
problem().temperatureAtPos(globalPos),
p_, pseudoFluidState);
Scalar v_n_ = 1. / FluidSystem::density(pseudoFluidState, nPhaseIdx);
problem().variables().dv_dp(globalIdxI) = (problem().variables().totalConcentration(globalIdxI, wCompIdx)
+ problem().variables().totalConcentration(globalIdxI, nCompIdx))
* ( v_n_ - 1./densityNWI) / -incp;
......@@ -1101,9 +1092,9 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
// initialize viscosities
problem().variables().viscosityWetting(globalIdx)
= FluidSystem::phaseViscosity(wPhaseIdx, temperature_, pressure[wPhaseIdx], fluidState);
= FluidSystem::viscosity(fluidState, wPhaseIdx);
problem().variables().viscosityNonwetting(globalIdx)
= FluidSystem::phaseViscosity(nPhaseIdx, temperature_, pressure[nPhaseIdx], fluidState);
= FluidSystem::viscosity(fluidState, nPhaseIdx);
// initialize mobilities
problem().variables().mobilityWetting(globalIdx) =
......@@ -1114,14 +1105,14 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
/ problem().variables().viscosityNonwetting(globalIdx);
// initialize mass fractions
problem().variables().wet_X1(globalIdx) = fluidState.massFrac(wPhaseIdx, wCompIdx);
problem().variables().nonwet_X1(globalIdx) = fluidState.massFrac(nPhaseIdx, wCompIdx);
problem().variables().wet_X1(globalIdx) = fluidState.massFraction(wPhaseIdx, wCompIdx);
problem().variables().nonwet_X1(globalIdx) = fluidState.massFraction(nPhaseIdx, wCompIdx);
// initialize densities
problem().variables().densityWetting(globalIdx)
= FluidSystem::phaseDensity(wPhaseIdx, temperature_, pressure[wPhaseIdx], fluidState);
= FluidSystem::density(fluidState, wPhaseIdx);
problem().variables().densityNonwetting(globalIdx)
= FluidSystem::phaseDensity(nPhaseIdx, temperature_, pressure[nPhaseIdx], fluidState);
= FluidSystem::density(fluidState, nPhaseIdx);
// determine volume mismatch between actual fluid volume and pore volume
Scalar sumConc = (problem().variables().totalConcentration(globalIdx, wCompIdx)
......@@ -1173,13 +1164,13 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
{
// Todo: only variables of present phase need to be updated
Scalar press1p = problem().variables().pressure()[globalIdx];
pseudoFluidState.update(Z1, press1p, problem().variables().saturation(globalIdx));
pseudoFluidState.update(Z1, press1p, problem().variables().saturation(globalIdx), problem().temperatureAtPos(globalPos));
// initialize viscosities
problem().variables().viscosityWetting(globalIdx)
= FluidSystem::phaseViscosity(wPhaseIdx, temperature_, press1p, pseudoFluidState);
= FluidSystem::viscosity(pseudoFluidState, wPhaseIdx);
problem().variables().viscosityNonwetting(globalIdx)
= FluidSystem::phaseViscosity(nPhaseIdx, temperature_, press1p, pseudoFluidState);
= FluidSystem::viscosity(pseudoFluidState, nPhaseIdx);
// initialize mobilities
problem().variables().mobilityWetting(globalIdx) =
......@@ -1190,14 +1181,14 @@ void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws()
/ problem().variables().viscosityNonwetting(globalIdx);
// initialize mass fractions
problem().variables().wet_X1(globalIdx) = pseudoFluidState.massFrac(wPhaseIdx, wCompIdx);
problem().variables().nonwet_X1(globalIdx) = pseudoFluidState.massFrac(nPhaseIdx, wCompIdx);
problem().variables().wet_X1(globalIdx) = pseudoFluidState.massFraction(wPhaseIdx, wCompIdx);
problem().variables().nonwet_X1(globalIdx) = pseudoFluidState.massFraction(nPhaseIdx, wCompIdx);
// initialize densities
problem().variables().densityWetting(globalIdx)
= FluidSystem::phaseDensity(wPhaseIdx, temperature_, press1p, pseudoFluidState);
= FluidSystem::density(pseudoFluidState, wPhaseIdx);
problem().variables().densityNonwetting(globalIdx)
= FluidSystem::phaseDensity(nPhaseIdx, temperature_, press1p, pseudoFluidState);
= FluidSystem::density(pseudoFluidState, nPhaseIdx);
// error term handling
Scalar sumConc = (problem().variables().totalConcentration(globalIdx, wCompIdx)
......
......@@ -365,16 +365,12 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
pressBound);
// determine fluid properties at the boundary
Scalar Xw1Bound = BCfluidState.massFrac(wPhaseIdx, wCompIdx);
Scalar Xn1Bound = BCfluidState.massFrac(nPhaseIdx, wCompIdx);
Scalar Xw1Bound = BCfluidState.massFraction(wPhaseIdx, wCompIdx);
Scalar Xn1Bound = BCfluidState.massFraction(nPhaseIdx, wCompIdx);
Scalar densityWBound = BCfluidState.density(wPhaseIdx);
Scalar densityNWBound = BCfluidState.density(nPhaseIdx);
Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx,
problem_.temperatureAtPos(globalPosFace),
pressBound[wPhaseIdx], BCfluidState);
Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx,
problem_.temperatureAtPos(globalPosFace),
pressBound[wPhaseIdx], BCfluidState);
Scalar viscosityWBound = FluidSystem::viscosity(BCfluidState, wPhaseIdx);
Scalar viscosityNWBound = FluidSystem::viscosity(BCfluidState, nPhaseIdx);
if(GET_PROP_VALUE(TypeTag, PTAG(EnableCapillarity)))
pcBound = BCfluidState.capillaryPressure();
......
......@@ -46,13 +46,18 @@ class PseudoOnePTwoCFluidState : public FluidState<typename GET_PROP_TYPE(TypeTa
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
public:
enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),
wPhaseIdx = 0,
nPhaseIdx = 1,};
numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents))};
enum {
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
wCompIdx = Indices::wPhaseIdx,
nCompIdx = Indices::nPhaseIdx,
};
public:
/*!
......@@ -70,31 +75,32 @@ public:
* \param press1p Pressure value for present phase \f$\mathrm{[Pa]}\f$
* \param satW Saturation of the wetting phase \f$\mathrm{[-]}\f$
*/
void update(const Scalar& Z1,const Scalar& press1p,const Scalar& satW)
void update(const Scalar& Z1,const Scalar& press1p,const Scalar& satW, const Scalar& temperature)
{
Sw_ = satW;
pressure1p_=press1p;
temperature_ = temperature;
if (satW == 1.)
{
massfracX1_[wPhaseIdx] = Z1;
massfracX1_[nPhaseIdx] = 0.;
massFractionWater_[wPhaseIdx] = Z1;
massFractionWater_[nPhaseIdx] = 0.;
molefracx1_[wPhaseIdx] = Z1 / FluidSystem::molarMass(0);
molefracx1_[wPhaseIdx] /= ( Z1 / FluidSystem::molarMass(0)
moleFractionWater_[wPhaseIdx] = Z1 / FluidSystem::molarMass(0);
moleFractionWater_[wPhaseIdx] /= ( Z1 / FluidSystem::molarMass(0)
+ (1.-Z1) / FluidSystem::molarMass(1));
molefracx1_[nPhaseIdx] = 0.;
moleFractionWater_[nPhaseIdx] = 0.;
}
else if (satW == 0.)
{
massfracX1_[wPhaseIdx] = 0.;
massfracX1_[nPhaseIdx] = Z1;
massFractionWater_[wPhaseIdx] = 0.;
massFractionWater_[nPhaseIdx] = Z1;
// interested in nComp => 1-X1
molefracx1_[nPhaseIdx] = ( Z1 / FluidSystem::molarMass(0) ); // = moles of compIdx
molefracx1_[nPhaseIdx] /= (Z1/ FluidSystem::molarMass(0)
moleFractionWater_[nPhaseIdx] = ( Z1 / FluidSystem::molarMass(0) ); // = moles of compIdx
moleFractionWater_[nPhaseIdx] /= (Z1/ FluidSystem::molarMass(0)
+ (1.-Z1) / FluidSystem::molarMass(1) ); // /= total moles in phase
molefracx1_[nPhaseIdx] = 0.;
moleFractionWater_[nPhaseIdx] = 0.;
}
else
Dune::dgrave << "Twophase conditions in single-phase flash! Saturation is " << satW << std::endl;
......@@ -118,7 +124,7 @@ public:
/*! \brief Returns the pressure of a fluid phase \f$\mathrm{[Pa]}\f$.
* \param phaseIdx Index of the phase
*/
Scalar phasePressure(int phaseIdx) const
Scalar pressure(int phaseIdx) const
{ return pressure1p_; }
......@@ -127,12 +133,12 @@ public:
* \param phaseIdx Index of the phase
* \param compIdx Index of the component
*/
Scalar massFrac(int phaseIdx, int compIdx) const
Scalar massFraction(int phaseIdx, int compIdx) const
{
if (compIdx == wPhaseIdx)
return massfracX1_[phaseIdx];
return massFractionWater_[phaseIdx];
else
return 1.-massfracX1_[phaseIdx];
return 1.-massFractionWater_[phaseIdx];
}
......@@ -141,20 +147,36 @@ public:
* \param phaseIdx Index of the phase
* \param compIdx Index of the component
*/
Scalar moleFrac(int phaseIdx, int compIdx) const
Scalar moleFraction(int phaseIdx, int compIdx) const
{
if (compIdx == wPhaseIdx)
return molefracx1_[phaseIdx];
return moleFractionWater_[phaseIdx];
else
return 1.-molefracx1_[phaseIdx];
return 1.-moleFractionWater_[phaseIdx];
}
Scalar averageMolarMass(int phaseIdx) const
{
return moleFractionWater_[phaseIdx]*FluidSystem::molarMass(wCompIdx)
+ (1.-moleFractionWater_[phaseIdx])*FluidSystem::molarMass(nCompIdx);
}
/*!
* \brief Returns the temperature of the fluids \f$\mathrm{[K]}\f$.
*
* Note that we assume thermodynamic equilibrium, so all fluids
* and the rock matrix exhibit the same temperature.
*/
Scalar temperature(int phaseIdx) const
{ return temperature_; };
//@}
public:
Scalar massfracX1_[numPhases];
Scalar molefracx1_[numPhases];
Scalar massFractionWater_[numPhases];
Scalar moleFractionWater_[numPhases];
Scalar Sw_;
Scalar pressure1p_;
Scalar temperature_;
};