Commit 58412316 authored by Beatrix Becker's avatar Beatrix Becker
Browse files

[2p2c sequential] calculate viscosity in flash, cleanup

parent a78e605e
...@@ -154,6 +154,9 @@ public: ...@@ -154,6 +154,9 @@ public:
fluidState.setMolarDensity(phase0Idx, FluidSystem::molarDensity(fluidState, phase0Idx)); fluidState.setMolarDensity(phase0Idx, FluidSystem::molarDensity(fluidState, phase0Idx));
fluidState.setMolarDensity(phase1Idx, FluidSystem::molarDensity(fluidState, phase1Idx)); fluidState.setMolarDensity(phase1Idx, FluidSystem::molarDensity(fluidState, phase1Idx));
fluidState.setViscosity(phase0Idx, FluidSystem::viscosity(fluidState, phase0Idx));
fluidState.setViscosity(phase1Idx, FluidSystem::viscosity(fluidState, phase1Idx));
Scalar sw = phaseMassFraction[phase0Idx] / fluidState.density(phase0Idx); Scalar sw = phaseMassFraction[phase0Idx] / fluidState.density(phase0Idx);
sw /= (phaseMassFraction[phase0Idx] / fluidState.density(phase0Idx) sw /= (phaseMassFraction[phase0Idx] / fluidState.density(phase0Idx)
+ phaseMassFraction[phase1Idx] / fluidState.density(phase1Idx)); + phaseMassFraction[phase1Idx] / fluidState.density(phase1Idx));
...@@ -185,10 +188,13 @@ public: ...@@ -185,10 +188,13 @@ public:
fluidState.setPresentPhaseIdx(presentPhaseIdx); fluidState.setPresentPhaseIdx(presentPhaseIdx);
fluidState.setMassFraction(presentPhaseIdx,comp0Idx, Z0); fluidState.setMassFraction(presentPhaseIdx,comp0Idx, Z0);
fluidState.setMassFraction(presentPhaseIdx,comp1Idx, 1. - Z0);
// transform mass to mole fractions // transform mass to mole fractions
fluidState.setMoleFraction(presentPhaseIdx, comp0Idx, Z0 / FluidSystem::molarMass(comp0Idx) fluidState.setMoleFraction(presentPhaseIdx, comp0Idx, Z0 / FluidSystem::molarMass(comp0Idx)
/ (Z0 / FluidSystem::molarMass(comp0Idx) + (1. - Z0) / FluidSystem::molarMass(comp1Idx))); / (Z0 / FluidSystem::molarMass(comp0Idx) + (1. - Z0) / FluidSystem::molarMass(comp1Idx)));
fluidState.setMoleFraction(presentPhaseIdx, comp1Idx, (1. - Z0) / FluidSystem::molarMass(comp1Idx)
/ (Z0 / FluidSystem::molarMass(comp0Idx) + (1. - Z0) / FluidSystem::molarMass(comp1Idx)));
fluidState.setAverageMolarMass(presentPhaseIdx, fluidState.setAverageMolarMass(presentPhaseIdx,
fluidState.massFraction(presentPhaseIdx, comp0Idx) * FluidSystem::molarMass(comp0Idx) fluidState.massFraction(presentPhaseIdx, comp0Idx) * FluidSystem::molarMass(comp0Idx)
...@@ -196,6 +202,8 @@ public: ...@@ -196,6 +202,8 @@ public:
fluidState.setDensity(presentPhaseIdx, FluidSystem::density(fluidState, presentPhaseIdx)); fluidState.setDensity(presentPhaseIdx, FluidSystem::density(fluidState, presentPhaseIdx));
fluidState.setMolarDensity(presentPhaseIdx, FluidSystem::molarDensity(fluidState, presentPhaseIdx)); fluidState.setMolarDensity(presentPhaseIdx, FluidSystem::molarDensity(fluidState, presentPhaseIdx));
fluidState.setViscosity(presentPhaseIdx, FluidSystem::viscosity(fluidState, presentPhaseIdx));
} }
//@} //@}
...@@ -254,6 +262,9 @@ public: ...@@ -254,6 +262,9 @@ public:
fluidState.setMolarDensity(phase0Idx, FluidSystem::molarDensity(fluidState, phase0Idx)); fluidState.setMolarDensity(phase0Idx, FluidSystem::molarDensity(fluidState, phase0Idx));
fluidState.setMolarDensity(phase1Idx, FluidSystem::molarDensity(fluidState, phase1Idx)); fluidState.setMolarDensity(phase1Idx, FluidSystem::molarDensity(fluidState, phase1Idx));
fluidState.setViscosity(phase0Idx, FluidSystem::viscosity(fluidState, phase0Idx));
fluidState.setViscosity(phase1Idx, FluidSystem::viscosity(fluidState, phase1Idx));
// set saturation // set saturation
fluidState.setSaturation(phase0Idx, saturation); fluidState.setSaturation(phase0Idx, saturation);
fluidState.setSaturation(phase1Idx, 1.0-saturation); fluidState.setSaturation(phase1Idx, 1.0-saturation);
......
...@@ -901,13 +901,13 @@ void FVPressure2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& en ...@@ -901,13 +901,13 @@ void FVPressure2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& en
* \param postTimeStep Flag indicating if we have just completed a time step * \param postTimeStep Flag indicating if we have just completed a time step
*/ */
template<class TypeTag> template<class TypeTag>
void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& elementI, bool postTimeStep) void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element, bool postTimeStep)
{ {
// get global coordinate of cell center // get global coordinate of cell center
GlobalPosition globalPos = elementI.geometry().center(); GlobalPosition globalPos = element.geometry().center();
// cell Index and cell data // cell Index and cell data
int eIdxGlobal = problem().variables().index(elementI); int eIdxGlobal = problem().variables().index(element);
CellData& cellData = problem().variables().cellData(eIdxGlobal); CellData& cellData = problem().variables().cellData(eIdxGlobal);
// acess the fluid state and prepare for manipulation // acess the fluid state and prepare for manipulation
...@@ -924,7 +924,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element ...@@ -924,7 +924,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
/ (cellData.massConcentration(wCompIdx) / (cellData.massConcentration(wCompIdx)
+ cellData.massConcentration(nCompIdx)); + cellData.massConcentration(nCompIdx));
// make shure only physical quantities enter flash calculation // make sure only physical quantities enter flash calculation
if(Z0 < 0. || Z0 > 1.) if(Z0 < 0. || Z0 > 1.)
{ {
Dune::dgrave << "Feed mass fraction unphysical: Z0 = " << Z0 Dune::dgrave << "Feed mass fraction unphysical: Z0 = " << Z0
...@@ -951,84 +951,60 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element ...@@ -951,84 +951,60 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
} }
} }
//determine phase pressures from primary pressure variable and pc of last TS PhaseVector pressure;
PhaseVector pressure(0.);
switch (pressureType)
{
case pw:
{
pressure[wPhaseIdx] = asImp_().pressure(eIdxGlobal);
pressure[nPhaseIdx] = asImp_().pressure(eIdxGlobal)
+ cellData.capillaryPressure();
break;
}
case pn:
{
pressure[wPhaseIdx] = asImp_().pressure(eIdxGlobal)
- cellData.capillaryPressure();
pressure[nPhaseIdx] = asImp_().pressure(eIdxGlobal);
break;
}
}
//complete fluid state
CompositionalFlash<Scalar, FluidSystem> flashSolver; CompositionalFlash<Scalar, FluidSystem> flashSolver;
flashSolver.concentrationFlash2p2c(fluidState, Z0, pressure, problem().spatialParams().porosity(elementI), temperature_); if(GET_PROP_VALUE(TypeTag, EnableCapillarity)) // iterate capillary pressure and saturation
// iterations part in case of enabled capillary pressure
Scalar pc(0.), oldPc(0.);
if(GET_PROP_VALUE(TypeTag, EnableCapillarity))
{ {
pc = MaterialLaw::pc(problem().spatialParams().materialLawParams(elementI), unsigned int maxiter = 6;
fluidState.saturation(wPhaseIdx)); Scalar pc = cellData.capillaryPressure(); // initial guess for pc from last TS
unsigned int maxiter = 5; // start iteration loop
int iterout = -1;
//start iteration loop
for (unsigned int iter = 0; iter < maxiter; iter++) for (unsigned int iter = 0; iter < maxiter; iter++)
{ {
//prepare pressures to enter flash calculation
switch (pressureType) switch (pressureType)
{ {
case pw: case pw:
{ {
// pressure[w] does not change, since it is primary variable pressure[wPhaseIdx] = asImp_().pressure(eIdxGlobal);
pressure[nPhaseIdx] = pressure[wPhaseIdx] + pc; pressure[nPhaseIdx] = asImp_().pressure(eIdxGlobal) + pc;
break; break;
} }
case pn: case pn:
{ {
// pressure[n] does not change, since it is primary variable pressure[wPhaseIdx] = asImp_().pressure(eIdxGlobal) - pc;
pressure[wPhaseIdx] = pressure[nPhaseIdx] - pc; pressure[nPhaseIdx] = asImp_().pressure(eIdxGlobal);
break; break;
} }
} }
//store old pc // complete fluid state
oldPc = pc;
//update with better pressures
flashSolver.concentrationFlash2p2c(fluidState, Z0, pressure, flashSolver.concentrationFlash2p2c(fluidState, Z0, pressure,
problem().spatialParams().porosity(elementI), problem().temperatureAtPos(globalPos)); problem().spatialParams().porosity(element), temperature_);
pc = MaterialLaw::pc(problem().spatialParams().materialLawParams(elementI),
fluidState.saturation(wPhaseIdx)); // calculate new pc
// TODO: get right criterion, do output for evaluation Scalar oldPc = pc;
//converge criterion pc = MaterialLaw::pc(problem().spatialParams().materialLawParams(element),
if (fabs(oldPc-pc)<10) fluidState.saturation(wPhaseIdx));
maxiter = 1;
iterout = iter; if (fabs(oldPc-pc)<10 && iter != 0)
break;
if (iter++ == maxiter)
Dune::dinfo << iter << "times iteration of pc was applied at Idx " << eIdxGlobal
<< ", pc delta still " << fabs(oldPc-pc) << std::endl;
} }
if (iterout != 0)
Dune::dinfo << iterout << "times iteration of pc was applied at Idx " << eIdxGlobal
<< ", pc delta still " << fabs(oldPc-pc) << std::endl;
} }
// initialize phase properties not stored in fluidstate else // capillary pressure neglected
cellData.setViscosity(wPhaseIdx, FluidSystem::viscosity(fluidState, wPhaseIdx)); {
cellData.setViscosity(nPhaseIdx, FluidSystem::viscosity(fluidState, nPhaseIdx)); pressure[wPhaseIdx] = pressure[nPhaseIdx] = asImp_().pressure()[eIdxGlobal];
flashSolver.concentrationFlash2p2c(fluidState, Z0, pressure,
problem().spatialParams().porosity(element), temperature_);
}
// initialize mobilities // initialize mobilities
cellData.setMobility(wPhaseIdx, MaterialLaw::krw(problem().spatialParams().materialLawParams(elementI), cellData.setMobility(wPhaseIdx, MaterialLaw::krw(problem().spatialParams().materialLawParams(element),
fluidState.saturation(wPhaseIdx)) fluidState.saturation(wPhaseIdx))
/ cellData.viscosity(wPhaseIdx)); / cellData.viscosity(wPhaseIdx));
cellData.setMobility(nPhaseIdx, MaterialLaw::krn(problem().spatialParams().materialLawParams(elementI), cellData.setMobility(nPhaseIdx, MaterialLaw::krn(problem().spatialParams().materialLawParams(element),
fluidState.saturation(wPhaseIdx)) fluidState.saturation(wPhaseIdx))
/ cellData.viscosity(nPhaseIdx)); / cellData.viscosity(nPhaseIdx));
...@@ -1043,7 +1019,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element ...@@ -1043,7 +1019,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
Scalar vol = massw / cellData.density(wPhaseIdx) + massn / cellData.density(nPhaseIdx); Scalar vol = massw / cellData.density(wPhaseIdx) + massn / cellData.density(nPhaseIdx);
if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(problem().timeManager().timeStepSize(), 0.0, 1.0e-30)) if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(problem().timeManager().timeStepSize(), 0.0, 1.0e-30))
{ {
cellData.volumeError()=(vol - problem().spatialParams().porosity(elementI)); cellData.volumeError()=(vol - problem().spatialParams().porosity(element));
using std::isnan; using std::isnan;
if (isnan(cellData.volumeError())) if (isnan(cellData.volumeError()))
...@@ -1052,7 +1028,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element ...@@ -1052,7 +1028,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
<< "volErr[" << eIdxGlobal << "] isnan: vol = " << vol << "volErr[" << eIdxGlobal << "] isnan: vol = " << vol
<< ", massw = " << massw << ", rho_l = " << cellData.density(wPhaseIdx) << ", massw = " << massw << ", rho_l = " << cellData.density(wPhaseIdx)
<< ", massn = " << massn << ", rho_g = " << cellData.density(nPhaseIdx) << ", massn = " << massn << ", rho_g = " << cellData.density(nPhaseIdx)
<< ", poro = " << problem().spatialParams().porosity(elementI) << ", poro = " << problem().spatialParams().porosity(element)
<< ", dt = " << problem().timeManager().timeStepSize()); << ", dt = " << problem().timeManager().timeStepSize());
} }
} }
......
...@@ -134,7 +134,7 @@ public: ...@@ -134,7 +134,7 @@ public:
Scalar dt_estimate = 0.; Scalar dt_estimate = 0.;
Dune::dinfo << "secant guess"<< std::endl; Dune::dinfo << "secant guess"<< std::endl;
problem_.transportModel().update(problem_.timeManager().time(), dt_estimate, updateEstimate_, false); problem_.transportModel().update(problem_.timeManager().time(), dt_estimate, updateEstimate_, false);
//last argument false in update() makes shure that this is estimate and no "real" transport step //last argument false in update() makes sure that this is estimate and no "real" transport step
// if we just started a new episode, the TS size of the update Estimate is a better // if we just started a new episode, the TS size of the update Estimate is a better
// estimate then the size of the last time step // estimate then the size of the last time step
...@@ -670,10 +670,6 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional) ...@@ -670,10 +670,6 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
problem_.transportModel().totalConcentration(wCompIdx,eIdxGlobal) = cellData.massConcentration(wCompIdx); problem_.transportModel().totalConcentration(wCompIdx,eIdxGlobal) = cellData.massConcentration(wCompIdx);
problem_.transportModel().totalConcentration(nCompIdx,eIdxGlobal) = cellData.massConcentration(nCompIdx); problem_.transportModel().totalConcentration(nCompIdx,eIdxGlobal) = cellData.massConcentration(nCompIdx);
// initialize phase properties not stored in fluidstate
cellData.setViscosity(wPhaseIdx, FluidSystem::viscosity(fluidState, wPhaseIdx));
cellData.setViscosity(nPhaseIdx, FluidSystem::viscosity(fluidState, nPhaseIdx));
// initialize mobilities // initialize mobilities
cellData.setMobility(wPhaseIdx, MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), cellData.setMobility(wPhaseIdx, MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
fluidState.saturation(wPhaseIdx)) fluidState.saturation(wPhaseIdx))
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment