Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Mixing energy #309

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Mixing energy #309

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
4c45158
changes to make black oil energy model less dependent not gennereting…
hnil Oct 17, 2018
1bc7c51
moved back to include enthapy
hnil Oct 17, 2018
f261e12
changed according to pullrequest comments
hnil Nov 6, 2018
df1d0a7
fixed error with template conversion
hnil Nov 7, 2018
cb99ef7
Merge branch 'mixing_energy' of github.com:hnil/opm-material into mix…
hnil Nov 7, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 7 additions & 1 deletion opm/material/fluidstates/BlackOilFluidState.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -321,7 +321,13 @@ class BlackOilFluidState
* exception!
*/
Scalar internalEnergy(unsigned phaseIdx OPM_UNUSED) const
{ return (*enthalpy_)[phaseIdx] - pressure(phaseIdx)/density(phaseIdx); }
{
//double entFac=0.0;
double entFac=1.0;
// this may be setting which for entFac =1 make this probably E100 compatible
// however this has to be consitent with call to enthalpy
return (*enthalpy_)[phaseIdx] - entFac*pressure(phaseIdx)/density(phaseIdx);
}

//////
// slow methods
Expand Down
102 changes: 75 additions & 27 deletions opm/material/fluidsystems/BlackOilFluidSystem.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -551,7 +551,7 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
if (enableDissolvedGas()) {
// miscible oil
const LhsEval& Rs = Opm::BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
const LhsEval& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
const LhsEval& bo = inverseFormationVolumeFactor(fluidState, phaseIdx, regionIdx);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Isn't the new version identical to the old one? if no, the Rs variable above should be unused?!

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do not think so, since since the member function do the bookkeeping with of finding if the fluid is saturated or not. It probably had been better to take all of this directly from the state, since this will ensure consistency with the flow problem which is essential for physical behavoir.. (Else the mass in the energy equation is not consitent with the mass int the system and tempreratur can be gennerated.)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How could I take it directly form the state. (I tried but broke some tests)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The interpolation between the saturated and the undersaturated interpolation should only affect the derivatives, i.e. if there are no convergence issues due to energy, it should not be necessary!?

How could I take it directly form the state. (I tried but broke some tests)

BlackOilFluidState has and invB() method. for generic fluid states this only makes limited sense, so it needs to be calculated. The best technical solution is probably to create something analogous to the getRs_(fluidState) function for inverse formation volume factors:

https://github.com/OPM/opm-material/blob/master/opm/material/fluidsystems/BlackOilFluidSystem.hpp#L52

But I guess that this is too much effort here...


return
bo*referenceDensity(oilPhaseIdx, regionIdx)
Expand All @@ -560,7 +560,7 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S

// immiscible oil
const LhsEval Rs(0.0);
const auto& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
const auto& bo = inverseFormationVolumeFactor(fluidState, phaseIdx, regionIdx);

return referenceDensity(phaseIdx, regionIdx)*bo;
}
Expand All @@ -569,16 +569,15 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
if (enableVaporizedOil()) {
// miscible gas
const LhsEval& Rv = Opm::BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);

const LhsEval& bg = inverseFormationVolumeFactor(fluidState, phaseIdx, regionIdx);
return
bg*referenceDensity(gasPhaseIdx, regionIdx)
+ Rv*bg*referenceDensity(oilPhaseIdx, regionIdx);
}

// immiscible gas
const LhsEval Rv(0.0);
const auto& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
const auto& bg = inverseFormationVolumeFactor(fluidState, phaseIdx, regionIdx);
return bg*referenceDensity(phaseIdx, regionIdx);
}

Expand Down Expand Up @@ -655,6 +654,74 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
}

template <class FluidState, class LhsEval = typename FluidState::Scalar>
static LhsEval energy(const FluidState& fluidState,
unsigned phaseIdx,
unsigned regionIdx)
{
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

that function should be named either internalEnergy() or enthalpy(). Also, internalEnergy() and enthalpy() need to be consistent, i.e. h_alpha = u_alpha + p_alpha/rho_alpha.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I followed your example and this energy is only used to calculate enthalpy and could probably be made private.

Copy link
Contributor

@andlaus andlaus Oct 18, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ah, okay. maybe this can even replace the enthalpy() and/or internalEnergy() functions directly?

assert(0 <= phaseIdx && phaseIdx <= numPhases);
assert(0 <= regionIdx && regionIdx <= numRegions());

//const auto& p = fluidState.pressure(phaseIdx);
//const auto& T = fluidState.temperature(phaseIdx);
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
switch (phaseIdx) {
case oilPhaseIdx: {
if (enableDissolvedGas()) {
// miscible oil
const LhsEval& Rs = Opm::BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
//const LhsEval& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
const LhsEval& bo = inverseFormationVolumeFactor(fluidState, phaseIdx, regionIdx);
//const LhsEval& bo = Opm::decay<LhsEval>(fluidState.invB(phaseIdx));
const LhsEval& Rv = Opm::BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
LhsEval energy =
bo*referenceDensity(oilPhaseIdx, regionIdx)
*oilPvt_->internalEnergy(regionIdx, T, p, Rs)
+ Rs*bo*referenceDensity(gasPhaseIdx, regionIdx)
*gasPvt_->internalEnergy(regionIdx, T, p, Rv);//this a to simple model
return energy/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
}

// immiscible oil
const LhsEval Rs(0.0);
//const LhsEval& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
//const LhsEval& bo = Opm::decay<LhsEval>(fluidState.invB(phaseIdx));
return oilPvt_->internalEnergy(regionIdx, T, p, Rs);
}

case gasPhaseIdx: {
if (enableVaporizedOil()) {
// miscible gas
const LhsEval& Rs = Opm::BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
const LhsEval& Rv = Opm::BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
const LhsEval& bg = inverseFormationVolumeFactor(fluidState, phaseIdx, regionIdx);
//const LhsEval& bg = Opm::decay<LhsEval>(fluidState.invB(phaseIdx));
LhsEval energy =
bg*referenceDensity(gasPhaseIdx, regionIdx)
*gasPvt_->internalEnergy(regionIdx, T, p, Rv)
+ Rv*bg*referenceDensity(oilPhaseIdx, regionIdx)
*oilPvt_->internalEnergy(regionIdx, T, p, Rs);// this is a to simple model
return energy/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
}

// immiscible gas
const LhsEval Rv(0.0);
// const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);

return gasPvt_->internalEnergy(regionIdx, T, p, Rv);

}

case waterPhaseIdx:
return waterPvt_->internalEnergy(regionIdx, T, p);

}

throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
}


/*!
* \brief Returns the formation volume factor \f$B_\alpha\f$ of an "undersaturated"
* fluid phase
Expand Down Expand Up @@ -963,28 +1030,9 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
assert(0 <= regionIdx && regionIdx <= numRegions());

const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));

switch (phaseIdx) {
case oilPhaseIdx:
return
oilPvt_->internalEnergy(regionIdx, T, p, Opm::BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
+ p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);

case gasPhaseIdx:
return
gasPvt_->internalEnergy(regionIdx, T, p, Opm::BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
+ p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);

case waterPhaseIdx:
return
waterPvt_->internalEnergy(regionIdx, T, p)
+ p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);

default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
}

throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
double entFac = 1.0;
return energy<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx)
+ entFac*p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
}

/*!
Expand Down