updating the thp value if VFP table is valid

through function isVFPActive().
This commit is contained in:
Kai Bao 2018-11-13 13:11:05 +01:00
parent 0d1a4b2d13
commit 32b00b61f8
2 changed files with 42 additions and 52 deletions

View File

@ -331,7 +331,7 @@ namespace Opm
template <class ValueType>
ValueType calculateBhpFromThp(const std::vector<ValueType>& rates, const int control_index) const;
double calculateThpFromBhp(const std::vector<double>& rates, const int control_index, const double bhp) const;
double calculateThpFromBhp(const std::vector<double>& rates, const double bhp) const;
// get the mobility for specific perforation
void getMobility(const Simulator& ebosSimulator,

View File

@ -1084,43 +1084,36 @@ namespace Opm
StandardWell<TypeTag>::
updateThp(WellState& well_state) const
{
// for the wells having a THP constaint, we should update their thp value
// If it is under THP control, it will be set to be the target value.
// TODO: a better standard is probably whether we have the table to calculate the THP value
// TODO: it is something we need to check the output to decide.
const WellControls* wc = well_controls_;
// TODO: we should only maintain one current control either from the well_state or from well_controls struct.
// Either one can be more favored depending on the final strategy for the initilzation of the well control
const int nwc = well_controls_get_num(wc);
// Looping over all controls until we find a THP constraint
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(wc, ctrl_index) == THP) {
// the current control
const int current = well_state.currentControls()[index_of_well_];
// if well under THP control at the moment
if (current == ctrl_index) {
const double thp_target = well_controls_iget_target(wc, current);
well_state.thp()[index_of_well_] = thp_target;
} else { // otherwise we calculate the thp from the bhp value
const Opm::PhaseUsage& pu = phaseUsage();
std::vector<double> rates(3, 0.0);
// When there is no vaild VFP table provided, we set the thp to be zero.
if (!this->isVFPActive()) {
well_state.thp()[index_of_well_] = 0.;
return;
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Water ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Oil ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Gas ] ];
}
// avaiable VFP table is provided, we should update the thp value
const double bhp = well_state.bhp()[index_of_well_];
// if the well is under THP control, we should use its target value
if (well_controls_get_current_type(well_controls_) == THP) {
well_state.thp()[index_of_well_] = well_controls_get_current_target(well_controls_);
} else {
// the well is under other control types, we calculate the thp based on bhp and rates
std::vector<double> rates(3, 0.0);
well_state.thp()[index_of_well_] = calculateThpFromBhp(rates, ctrl_index, bhp);
}
break;
const Opm::PhaseUsage& pu = phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Water ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Oil ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Gas ] ];
}
// TODO: we should have some guard here to avoid some NAN or INF ratios
const double bhp = well_state.bhp()[index_of_well_];
well_state.thp()[index_of_well_] = calculateThpFromBhp(rates, bhp);
}
}
@ -2096,7 +2089,6 @@ namespace Opm
double
StandardWell<TypeTag>::
calculateThpFromBhp(const std::vector<double>& rates,
const int control_index,
const double bhp) const
{
assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
@ -2105,32 +2097,30 @@ namespace Opm
const double liquid = rates[Oil];
const double vapour = rates[Gas];
const int vfp = well_controls_iget_vfp(well_controls_, control_index);
const double& alq = well_controls_iget_alq(well_controls_, control_index);
// pick the density in the top layer
const double rho = perf_densities_[0];
double thp = 0.0;
if (well_type_ == INJECTOR) {
const double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth();
const int table_id = well_ecl_->getInjectionProperties(current_step_).VFPTableNumber;
const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
thp = vfp_properties_->getInj()->thp(vfp, aqua, liquid, vapour, bhp + dp);
}
else if (well_type_ == PRODUCER) {
const double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth();
thp = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp);
}
else if (well_type_ == PRODUCER) {
const int table_id = well_ecl_->getProductionProperties(current_step_).VFPTableNumber;
const double alq = well_ecl_->getProductionProperties(current_step_).ALQValue;
const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
thp = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq);
}
else {
OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
}
thp = vfp_properties_->getProd()->thp(vfp, aqua, liquid, vapour, bhp + dp, alq);
}
else {
OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
}
return thp;
return thp;
}