Commit 6dd98213 authored by Andreas Lauser's avatar Andreas Lauser
Browse files

peng-robinson test: make it work

this test can now be used to find the parameters for the black-oil model
(with slight modifications)

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7588 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 744be5e8
......@@ -184,10 +184,12 @@ public:
catch (Dune::FMatrixError e)
{
/*
printFluidState_(fluidState);
std::cout << "error: " << e << "\n";
std::cout << "b: " << b << "\n";
std::cout << "J: " << J << "\n";
*/
throw Dumux::NumericalProblem(e.what());
}
Valgrind::CheckDefined(deltaX);
......
......@@ -49,7 +49,6 @@ namespace Dumux
{
/*!
*
* \brief Implements the Peng-Robinson equation of state for liquids
* and gases.
*
......@@ -202,18 +201,22 @@ public:
else {
// find the extrema (if they are present)
Scalar Vmin, Vmax, pmin, pmax;
findExtrema_(Vmin, Vmax,
pmin, pmax,
a, b, T);
if (isGasPhase)
Vm = std::max(Vmax, VmCubic);
if (findExtrema_(Vmin, Vmax,
pmin, pmax,
a, b, T))
{
if (isGasPhase)
Vm = std::max(Vmax, VmCubic);
else
Vm = std::min(Vmin, VmCubic);
}
else
Vm = std::min(Vmin, VmCubic);
Vm = VmCubic;
}
}
Valgrind::CheckDefined(Vm);
assert(std::isfinite(Vm) && Vm > 0);
return Vm;
}
......@@ -418,7 +421,7 @@ protected:
// invert resulting cubic polynomial analytically
Scalar allV[4];
allV[0] = V;
int numSol = 1 + Dumux::invertCubicPolynomial(&allV[1], b1, b2, b3, b4);
int numSol = 1 + Dumux::invertCubicPolynomial(allV + 1, b1, b2, b3, b4);
// sort all roots of the derivative
std::sort(allV + 0, allV + numSol);
......
......@@ -93,6 +93,9 @@ public:
*/
void updatePure(Scalar temperature, Scalar pressure)
{
Valgrind::CheckDefined(temperature);
Valgrind::CheckDefined(pressure);
// Calculate the Peng-Robinson parameters of the pure
// components
//
......@@ -189,8 +192,7 @@ public:
*/
template <class FluidState>
void updateSingleMoleFraction(const FluidState &fs,
int compIdx,
Scalar delta)
int compIdx)
{
updateMix(fs);
}
......
......@@ -323,6 +323,7 @@ public:
int phaseIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
return fluidState.averageMolarMass(phaseIdx)/paramCache.molarVolume(phaseIdx);
}
......@@ -462,15 +463,13 @@ private:
case H2OIdx: return H2O::vaporPressure(temperature);
// the values of the Henry constant for the solutes have
// been computed using the Peng-Robinson equation of state
// (-> slope of the component's fugacity function at
// almost 100% water content)
// are faked so far...
case C1Idx: return 5.57601e+09;
case C3Idx: return 1.89465e+10;
case C6Idx: return 5.58969e+12;
case C10Idx: return 4.31947e+17;
case C15Idx: return 4.27283e+28;
case C20Idx: return 3.39438e+36;
case C3Idx: return 1e10;
case C6Idx: return 1e10;
case C10Idx: return 1e10;
case C15Idx: return 1e10;
case C20Idx: return 1e10;
default: DUNE_THROW(Dune::InvalidStateException, "Unknown component index " << compIdx);
}
};
......
......@@ -107,19 +107,9 @@ public:
int compIdx)
{
if (phaseIdx == oPhaseIdx)
oilPhaseParams_.updateSingleMoleFraction(fs,
compIdx,
fs.moleFraction(phaseIdx, compIdx)
- moleFrac_[phaseIdx][compIdx]);
oilPhaseParams_.updateSingleMoleFraction(fs, compIdx);
else if (phaseIdx == gPhaseIdx)
gasPhaseParams_.updateSingleMoleFraction(fs,
compIdx,
fs.moleFraction(phaseIdx, compIdx)
- moleFrac_[phaseIdx][compIdx]);
// update the mole fraction which the parameters are
// calculated for
moleFrac_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx);
gasPhaseParams_.updateSingleMoleFraction(fs, compIdx);
// update the phase's molar volume
updateMolarVolume_(fs, phaseIdx);
......@@ -226,20 +216,16 @@ public:
{
updatePure_(fs, phaseIdx);
updateMix_(fs, phaseIdx);
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
moleFrac_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx);
VmUpToDate_[phaseIdx] = false;
}
else if (!(exceptQuantities & ParentType::Composition))
{
updateMix_(fs, phaseIdx);
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
moleFrac_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx);
VmUpToDate_[phaseIdx] = false;
}
else if (!(exceptQuantities & ParentType::Pressure))
else if (!(exceptQuantities & ParentType::Pressure)) {
VmUpToDate_[phaseIdx] = false;
}
}
protected:
......@@ -340,7 +326,6 @@ protected:
bool VmUpToDate_[numPhases];
Scalar Vm_[numPhases];
Scalar moleFrac_[numPhases][numComponents];
OilPhaseParams oilPhaseParams_;
GasPhaseParams gasPhaseParams_;
......
......@@ -88,7 +88,7 @@ int main(int argc, char** argv)
fluidState.setPressure(oPhaseIdx, 4000 * 6894.7573); // 4000 PSI
// oil saturation
fluidState.setPressure(oPhaseIdx, 1.0);
fluidState.setSaturation(oPhaseIdx, 1.0);
// composition: SPE-5 reservoir oil
fluidState.setMoleFraction(oPhaseIdx, H2OIdx, 0.0);
......@@ -109,7 +109,7 @@ int main(int argc, char** argv)
////////////
ComponentVector molarities;
for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
molarities[compIdx] = fluidState.molarity(oPhaseIdx, compIdx);
molarities[compIdx] = fluidState.saturation(oPhaseIdx)*fluidState.molarity(oPhaseIdx, compIdx);
////////////
// Gradually increase the volume for and calculate the gas
......@@ -118,14 +118,15 @@ int main(int argc, char** argv)
////////////
FluidState flashFluidState;
flashFluidState.assign(fluidState);
Flash::guessInitial(flashFluidState, paramCache, molarities);
Flash::solve<MaterialLaw>(flashFluidState, paramCache, matParams, molarities);
std::cout << "p[Pa] rho_o[kg/m^3] rho_g[kg/m^3]\n";
std::cout << "alpha[-] p[Pa] S_g[-] rho_o[kg/m^3] rho_g[kg/m^3]\n";
int n = 1000;
for (int i = 0; i < n; ++i) {
Scalar minAlpha = 1.0;
Scalar maxAlpha = 3.0;
Scalar minAlpha = 0.98;
Scalar maxAlpha = 5.0;
// ratio between the original and the current volume
Scalar alpha = minAlpha + (maxAlpha - minAlpha)*i/(n - 1);
......@@ -147,7 +148,9 @@ int main(int argc, char** argv)
Scalar rhoGSurface = FluidSystem::density(surfaceFluidState, paramCache, gPhaseIdx);
Scalar rhoOSurface = FluidSystem::density(surfaceFluidState, paramCache, gPhaseIdx);
*/
std::cout << flashFluidState.pressure(oPhaseIdx) << " "
std::cout << alpha << " "
<< flashFluidState.pressure(oPhaseIdx) << " "
<< flashFluidState.saturation(gPhaseIdx) << " "
<< flashFluidState.density(oPhaseIdx) << " "
<< flashFluidState.density(gPhaseIdx) << " "
<< "\n";
......
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