ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim > Class Template Reference

Detailed Description

template<typename ShapeFunction, int GlobalDim>
class ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >

Definition at line 30 of file WellboreSimulatorFEM.h.

#include <WellboreSimulatorFEM.h>

Inheritance diagram for ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >:

Collaboration diagram for ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >:

Public Member Functions
	WellboreSimulatorFEM (MeshLib::Element const &element, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, WellboreSimulatorProcessData const &process_data)
void	assemble (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
Eigen::Map< const Eigen::RowVectorXd >	getShapeMatrix (const unsigned integration_point) const override
	Provides the shape matrix at the given integration point.
void	computeSecondaryVariableConcrete (double const, double const, Eigen::VectorXd const &, Eigen::VectorXd const &) override
void	postTimestepConcrete (Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override
Public Member Functions inherited from ProcessLib::WellboreSimulator::WellboreSimulatorLocalAssemblerInterface
	WellboreSimulatorLocalAssemblerInterface ()=default
Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual	~LocalAssemblerInterface ()=default
virtual void	setInitialConditions (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, int const process_id)
virtual void	initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
virtual void	preAssemble (double const, double const, std::vector< double > const &)
virtual void	assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
virtual void	assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
virtual void	assembleWithJacobianForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
virtual void	computeSecondaryVariable (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
virtual void	preTimestep (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
virtual void	postTimestep (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
void	postNonLinearSolver (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
virtual Eigen::Vector3d	getFlux (MathLib::Point3d const &, double const, std::vector< double > const &) const
virtual Eigen::Vector3d	getFlux (MathLib::Point3d const &, double const, std::vector< std::vector< double > > const &) const
	Fits to staggered scheme.
virtual int	getNumberOfVectorElementsForDeformation () const
Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual	~ExtrapolatableElement ()=default

Protected Types
using	IpData

Protected Member Functions
virtual std::vector< double > const &	getIntPtVaporMassFlowRate (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
virtual std::vector< double > const &	getIntPtLiquidMassFlowRate (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
virtual std::vector< double > const &	getIntPtTemperature (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
virtual std::vector< double > const &	getIntPtDryness (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
virtual std::vector< double > const &	getIntPtVaporVolumeFraction (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override

Protected Attributes
MeshLib::Element const &	_element
NumLib::GenericIntegrationMethod const &	_integration_method
bool const	_is_axially_symmetric
WellboreSimulatorProcessData const &	_process_data
Eigen::Vector3d	_element_direction
std::vector< IpData, Eigen::aligned_allocator< IpData > >	_ip_data

Static Protected Attributes
static const int	pressure_index = 0
static const int	velocity_index = ShapeFunction::NPOINTS
static const int	enthalpy_index = 2 * ShapeFunction::NPOINTS
static const int	pressure_size = ShapeFunction::NPOINTS
static const int	velocity_size = ShapeFunction::NPOINTS
static const int	enthalpy_size = ShapeFunction::NPOINTS

Private Types
using	ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>
using	ShapeMatrices = typename ShapeMatricesType::ShapeMatrices
using	LocalMatrixType
using	LocalVectorType
using	NodalVectorType = typename ShapeMatricesType::NodalVectorType
using	NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType
using	GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType
using	GlobalDimNodalMatrixType
using	GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType

Member Typedef Documentation

GlobalDimMatrixType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType

private

Definition at line 48 of file WellboreSimulatorFEM.h.

GlobalDimNodalMatrixType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::GlobalDimNodalMatrixType

private

Initial value:

 
        typename ShapeMatricesType::GlobalDimNodalMatrixType

Definition at line 46 of file WellboreSimulatorFEM.h.

GlobalDimVectorType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType

private

Definition at line 45 of file WellboreSimulatorFEM.h.

IpData

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::IpData

protected

Initial value:

 
        IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>

Definition at line 227 of file WellboreSimulatorFEM.h.

LocalMatrixType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::LocalMatrixType

private

Initial value:

 typename ShapeMatricesType::template MatrixType<
        NUM_NODAL_DOF * ShapeFunction::NPOINTS,
        NUM_NODAL_DOF * ShapeFunction::NPOINTS>

Definition at line 35 of file WellboreSimulatorFEM.h.

LocalVectorType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::LocalVectorType

private

Initial value:

 
        typename ShapeMatricesType::template VectorType<NUM_NODAL_DOF *
                                                        ShapeFunction::NPOINTS>

Definition at line 38 of file WellboreSimulatorFEM.h.

NodalRowVectorType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType

private

Definition at line 43 of file WellboreSimulatorFEM.h.

NodalVectorType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::NodalVectorType = typename ShapeMatricesType::NodalVectorType

private

Definition at line 42 of file WellboreSimulatorFEM.h.

ShapeMatrices

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::ShapeMatrices = typename ShapeMatricesType::ShapeMatrices

private

Definition at line 33 of file WellboreSimulatorFEM.h.

ShapeMatricesType

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>

private

Definition at line 32 of file WellboreSimulatorFEM.h.

Constructor & Destructor Documentation

WellboreSimulatorFEM()

template<typename ShapeFunction, int GlobalDim>

ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::WellboreSimulatorFEM ( MeshLib::Element const & element, std::size_t const , NumLib::GenericIntegrationMethod const & integration_method, bool const is_axially_symmetric, WellboreSimulatorProcessData const & process_data )

inline

Definition at line 51 of file WellboreSimulatorFEM.h.

        : _element(element),
          _integration_method(integration_method),
          _is_axially_symmetric(is_axially_symmetric),
          _process_data(process_data)
    {
        // calculate the element direction vector
        auto const& p0 = element.getNode(0)->asEigenVector3d();
        auto const& p1 = element.getNode(1)->asEigenVector3d();
 
        _element_direction = (p1 - p0).normalized();
 
        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();
        _ip_data.reserve(n_integration_points);
 
        auto const shape_matrices =
            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
                                      GlobalDim>(element, is_axially_symmetric,
                                                 _integration_method);
 
        ParameterLib::SpatialPosition pos;
        pos.setElementID(_element.getID());
        auto const& medium =
            *_process_data.media_map.getMedium(_element.getID());
        auto const& liquid_phase =
            medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
 
        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            _ip_data.emplace_back(
                shape_matrices[ip].N, shape_matrices[ip].dNdx,
                _integration_method.getWeightedPoint(ip).getWeight() *
                    shape_matrices[ip].integralMeasure *
                    shape_matrices[ip].detJ);
 
            MaterialPropertyLib::VariableArray vars;
 
            NodalVectorType ref_p =
                _process_data.well_ref_pressure.getNodalValuesOnElement(
                    _element, 0);
            NodalVectorType ref_h =
                _process_data.well_ref_enthalpy.getNodalValuesOnElement(
                    _element, 0);
 
            vars.liquid_phase_pressure = _ip_data[ip].N.dot(ref_p);
            vars.enthalpy = _ip_data[ip].N.dot(ref_h);
 
            //  .initialValue
            _ip_data[ip].temperature =
                liquid_phase
                    .property(MaterialPropertyLib::PropertyType::temperature)
                    .template value<double>(vars, pos, 0, 0);
            vars.temperature = _ip_data[ip].temperature;
            _ip_data[ip].mix_density =
                liquid_phase
                    .property(MaterialPropertyLib::PropertyType::density)
                    .template value<double>(vars, pos, 0, 0);
            _ip_data[ip].dryness = 0;
            _ip_data[ip].vapor_volume_fraction = 0;
            _ip_data[ip].vapor_mass_flow_rate = 0;
            _ip_data[ip].liquid_mass_flow_rate = 0;
            _ip_data[ip].pushBackState();
        }
    }

References _element, _element_direction, _integration_method, _ip_data, _is_axially_symmetric, _process_data, MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::enthalpy, MeshLib::Element::getNode(), NumLib::initShapeMatrices(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, ParameterLib::SpatialPosition::setElementID(), MaterialPropertyLib::temperature, and MaterialPropertyLib::VariableArray::temperature.

Member Function Documentation

assemble()

template<typename ShapeFunction, int GlobalDim>

void ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::assemble ( double const t, double const dt, std::vector< double > const & local_x, std::vector< double > const & local_x_prev, std::vector< double > & local_M_data, std::vector< double > & local_K_data, std::vector< double > & local_b_data )

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 31 of file WellboreSimulatorFEM-impl.h.

{
    auto const local_matrix_size = local_x.size();
 
    assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
 
    auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
        local_M_data, local_matrix_size, local_matrix_size);
    auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
        local_K_data, local_matrix_size, local_matrix_size);
    auto local_b = MathLib::createZeroedVector<LocalVectorType>(
        local_b_data, local_matrix_size);
 
    // Get block matrices
    auto Mvv = local_M.template block<velocity_size, velocity_size>(
        velocity_index, velocity_index);
 
    auto Mhp = local_M.template block<enthalpy_size, pressure_size>(
        enthalpy_index, pressure_index);
    auto Mhh = local_M.template block<enthalpy_size, enthalpy_size>(
        enthalpy_index, enthalpy_index);
 
    auto Kpv = local_K.template block<pressure_size, velocity_size>(
        pressure_index, velocity_index);
 
    auto Kvp = local_K.template block<velocity_size, pressure_size>(
        velocity_index, pressure_index);
    auto Kvv = local_K.template block<velocity_size, velocity_size>(
        velocity_index, velocity_index);
 
    auto Khh = local_K.template block<enthalpy_size, enthalpy_size>(
        enthalpy_index, enthalpy_index);
 
    auto Bp = local_b.template segment<pressure_size>(pressure_index);
    auto Bv = local_b.template segment<velocity_size>(velocity_index);
    auto Bh = local_b.template segment<enthalpy_size>(enthalpy_index);
 
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
 
    ParameterLib::SpatialPosition pos;
    pos.setElementID(_element.getID());
 
    auto const& b = _process_data.specific_body_force;
 
    MaterialPropertyLib::VariableArray vars;
 
    // Get material properties
    auto const& medium = *_process_data.media_map.getMedium(_element.getID());
    auto const& liquid_phase =
        medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
    auto const& gas_phase = medium.phase(MaterialPropertyLib::PhaseName::Gas);
 
    // Get wellbore parameters
    // casing thickness
    auto const t_ca = _process_data.wellbore.casing_thickness(t, pos)[0];
    // wellbore radius
    auto const r_w = _process_data.wellbore.diameter(t, pos)[0] / 2;
 
    // pipe thickness
    auto const t_p = _process_data.wellbore.pipe_thickness(t, pos)[0];
 
    // roughness of the wellbore
    auto const xi = _process_data.wellbore.roughness(t, pos)[0];
    // pipe outer radius
    auto const r_o = r_w - t_ca;
    // pipe inner radius
    auto const r_i = r_o - t_p;
 
    // get reservoir properties
    NodalVectorType T_r =
        _process_data.reservoir_properties.temperature.getNodalValuesOnElement(
            _element, t);
    NodalVectorType p_r =
        _process_data.reservoir_properties.pressure.getNodalValuesOnElement(
            _element, t);
    NodalVectorType PI =
        _process_data.productivity_index.getNodalValuesOnElement(_element, t);
    auto const k_r =
        _process_data.reservoir_properties.thermal_conductivity(t, pos)[0];
    auto const rho_r = _process_data.reservoir_properties.density(t, pos)[0];
    auto const c_r =
        _process_data.reservoir_properties.specific_heat_capacity(t, pos)[0];
 
    for (unsigned ip(0); ip < n_integration_points; ip++)
    {
        auto& ip_data = _ip_data[ip];
        auto const& N = ip_data.N;
        auto const& dNdx = ip_data.dNdx;
        auto const& w = ip_data.integration_weight;
        auto& mix_density = ip_data.mix_density;
        auto& temperature = ip_data.temperature;
        auto& steam_mass_frac = ip_data.dryness;
        auto& vapor_volume_frac = ip_data.vapor_volume_fraction;
        auto& vapor_mass_flowrate = ip_data.vapor_mass_flow_rate;
        auto& liquid_mass_flowrate = ip_data.liquid_mass_flow_rate;
 
        pos = {
            std::nullopt, _element.getID(),
            MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
                                                            ShapeMatricesType>(
                _element, N))};
 
        double p_int_pt = 0.0;
        double v_int_pt = 0.0;
        double h_int_pt = 0.0;
 
        NumLib::shapeFunctionInterpolate(local_x, N, p_int_pt, v_int_pt,
                                         h_int_pt);
 
        double p_prev_int_pt = 0.0;
        double v_prev_int_pt = 0.0;
        double h_prev_int_pt = 0.0;
 
        NumLib::shapeFunctionInterpolate(local_x_prev, N, p_prev_int_pt,
                                         v_prev_int_pt, h_prev_int_pt);
 
        double vdot_int_pt = (v_int_pt - v_prev_int_pt) / dt;
 
        // calculate fluid properties
 
        const double pi = std::numbers::pi;
 
        vars.liquid_phase_pressure = p_int_pt;
        vars.enthalpy = h_int_pt;
 
        // Above the critical pressure the region 4 saturation line ends, so
        // there is no two-phase state to describe and the saturation
        // properties are not evaluated at all: they would be extrapolated,
        // and the closure they feed has no admissible void fraction there.
        // Such a section is compressed liquid, which the region 1 properties
        // of the liquid phase describe, and it is solved as one. Below the
        // lower bound of the saturation line there is no such fall-back - the
        // water is vapour, which this process has no properties for - so the
        // range check of the saturation properties aborts the assembly as
        // before.
        double dryness = 0.;
        double T_int_pt = 0.;
        double liquid_water_density = 0.;
        double vapour_water_density = 0.;
        double alpha = 0.;
        std::optional<MaterialPropertyLib::DriftFluxState> drift_flux_state;
 
        if (p_int_pt >
            MaterialLib::PhysicalConstant::CriticalPoint::PressureWater)
        {
            T_int_pt =
                liquid_phase
                    .property(MaterialPropertyLib::PropertyType::temperature)
                    .template value<double>(vars, pos, t, dt);
            vars.temperature = T_int_pt;
 
            MaterialPropertyLib::IAPWSIF97Region1::checkStateInRange(
                p_int_pt, T_int_pt,
                "the compressed liquid state of the WellboreSimulator "
                "process");
 
            liquid_water_density =
                liquid_phase
                    .property(MaterialPropertyLib::PropertyType::density)
                    .template value<double>(vars, pos, t, dt);
        }
        else
        {
            liquid_water_density =
                liquid_phase
                    .property(
                        MaterialPropertyLib::PropertyType::saturation_density)
                    .template value<double>(vars, pos, t, dt);
            vapour_water_density =
                gas_phase
                    .property(
                        MaterialPropertyLib::PropertyType::saturation_density)
                    .template value<double>(vars, pos, t, dt);
 
            double const h_sat_liq_w =
                liquid_phase
                    .property(
                        MaterialPropertyLib::PropertyType::saturation_enthalpy)
                    .template value<double>(vars, pos, t, dt);
            double const h_sat_vap_w =
                gas_phase
                    .property(
                        MaterialPropertyLib::PropertyType::saturation_enthalpy)
                    .template value<double>(vars, pos, t, dt);
 
            dryness = MaterialPropertyLib::steamDryness(h_int_pt, h_sat_liq_w,
                                                        h_sat_vap_w);
 
            T_int_pt =
                (dryness == 0)
                    ? liquid_phase
                          .property(
                              MaterialPropertyLib::PropertyType::temperature)
                          .template value<double>(vars, pos, t, dt)
                    : gas_phase
                          .property(MaterialPropertyLib::PropertyType::
                                        saturation_temperature)
                          .template value<double>(vars, pos, t, dt);
            vars.temperature = T_int_pt;
 
            // For the calculation of the void fraction of vapour,
            // see Rohuani, Z., and E. Axelsson. "Calculation of volume void
            // fraction in a subcooled and quality region." International
            // Journal of Heat and Mass Transfer 17 (1970): 383-393.
 
            // The drift is aligned with the mixture flow so that the closure
            // below and the slip momentum term further down are consistent,
            // see MaterialPropertyLib::alignedDriftFluxVelocity().
            drift_flux_state = MaterialPropertyLib::driftFluxState(
                dryness, T_int_pt, vapour_water_density, liquid_water_density,
                v_int_pt);
 
            // solving void fraction of vapour: Rouhani-Axelsson
            auto const alpha_solution =
                MaterialPropertyLib::computeVapourVoidFraction(
                    *drift_flux_state);
 
            if (!alpha_solution)
            {
                throw NumLib::AssemblyException(fmt::format(
                    "The drift-flux closure of the WellboreSimulator process "
                    "has no admissible vapour void fraction in element {:d}, "
                    "integration point {:d}: pressure {:g} Pa, mixture "
                    "velocity {:g} m/s, specific enthalpy {:g} J/kg, "
                    "temperature {:g} K, {}",
                    _element.getID(), ip, p_int_pt, v_int_pt, h_int_pt,
                    T_int_pt,
                    MaterialPropertyLib::voidFractionClosureDiagnostics(
                        *drift_flux_state)));
            }
 
            alpha = *alpha_solution;
 
            if (alpha == 0)
            {
                liquid_water_density =
                    liquid_phase
                        .property(MaterialPropertyLib::PropertyType::density)
                        .template value<double>(vars, pos, t, dt);
            }
        }
 
        steam_mass_frac = dryness;
        temperature = T_int_pt;
        vapor_volume_frac = alpha;
 
        mix_density =
            vapour_water_density * alpha + liquid_water_density * (1 - alpha);
 
        auto& mix_density_prev = ip_data.mix_density_prev;
        vars.density = mix_density;
 
        auto const rho_dot = (mix_density - mix_density_prev) / dt;
 
        double const liquid_water_velocity_act =
            (alpha == 0)   ? v_int_pt
            : (alpha == 1) ? 0
                           : (1 - dryness) * mix_density * v_int_pt /
                                 (1 - alpha) / liquid_water_density;
        double const vapor_water_velocity_act =
            (alpha == 0) ? 0
                         : dryness * mix_density * v_int_pt /
                               (alpha * vapour_water_density);
 
        vapor_mass_flowrate = vapor_water_velocity_act * vapour_water_density *
                              pi * r_i * r_i * alpha;
 
        liquid_mass_flowrate = liquid_water_velocity_act *
                               liquid_water_density * pi * r_i * r_i *
                               (1 - alpha);
 
        double const gamma = drift_flux_state
                                 ? MaterialPropertyLib::mixtureSlipParameter(
                                       alpha, *drift_flux_state)
                                 : 0.;
 
        double const miu =
            liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity)
                .template value<double>(vars, pos, t, dt);
        double const Re = mix_density * v_int_pt * 2 * r_i / miu;
 
        // Wall friction coefficient,
        // Musivand Arzanfudi, Mehdi, and Rafid Al‐Khoury. "A compressible
        // two‐fluid multiphase model for CO2 leakage through a wellbore."
        // International Journal for Numerical Methods in Fluids 77.8 (2015):
        // 477-507.
        double f = 0.0;
        if (Re > 10 && Re <= 2400)
        {
            f = 16 / Re;
        }
        else if (Re > 2400)
        {
            f = std::pow(std::log(xi / 3.7 / r_i) -
                             5.02 / Re * std::log(xi / 3.7 / r_i + 13 / Re),
                         -2) /
                16;
        }
 
        double Q_hx = 0;
        double const T_r_int_pt = N.dot(T_r);
        // conductive heat exchange between wellbore and formation
        if (_process_data.has_heat_exchange_with_formation)
        {
            // See Zhang, Pan, Pruess, Finsterle (2011). A time-convolution
            // approach for modeling heat exchange between a wellbore and
            // surrounding formation. Geothermics 40, 261-266.
            const double alpha_r = k_r / rho_r / c_r;
            const double t_d = alpha_r * t / (r_i * r_i);
 
            double beta;
            if (t_d < 2.8)
            {
                beta = 1 / std::sqrt(pi * t_d) + 0.5 -
                       0.25 * std::sqrt(t_d / pi) + 0.125 * t_d;
            }
            else
            {
                beta = 2 * (1 / (std::log(4 * t_d) - 2 * 0.57722) -
                            0.57722 /
                                std::pow((std::log(4 * t_d) - 2 * 0.57722), 2));
            }
 
            const double P_c = 2 * pi * r_i;
            Q_hx = P_c * k_r * (T_r_int_pt - T_int_pt) / r_i * beta;
        }
 
        // mass exchange with reservoir
        double const p_r_int_pt = N.dot(p_r);
        double const PI_int_pt = N.dot(PI);
        double Q_mx = PI_int_pt * (p_int_pt - p_r_int_pt);
 
        // advective momentum and energy exchange with reservoir due to the mass
        // exchange
        double Q_mom = 0;
        double Q_ene = 0;
        if (Q_mx != 0)
        {
            Q_mom = Q_mx * v_int_pt;
            // Only single-phase liquid condition is considered now
            // TODO: update the two-phase flowing enthalpy from the feed zone.
            vars.liquid_phase_pressure = p_r_int_pt;
            vars.temperature = T_r_int_pt;
            double const h_fres =
                liquid_phase
                    .property(MaterialPropertyLib::PropertyType::enthalpy)
                    .template value<double>(vars, pos, t, dt);
            Q_ene = Q_mx * h_fres;
        }
 
        // M matrix assembly
        Mvv.noalias() += w * N.transpose() * mix_density * N;
 
        Mhp.noalias() += -w * N.transpose() * N;
        Mhh.noalias() += w * N.transpose() * mix_density * N;
 
        // K matrix assembly
        Kpv.noalias() += w * dNdx.transpose() * N * mix_density;
 
        Kvp.noalias() += w * N.transpose() * dNdx;
        Kvv.noalias() += w * N.transpose() * rho_dot * N;
 
        Khh.noalias() += w * N.transpose() * mix_density * v_int_pt * dNdx;
 
        // b matrix assembly
        Bp.noalias() += w * N.transpose() * rho_dot + w * N.transpose() * Q_mx;
 
        Bv.noalias() +=
            w * dNdx.transpose() * mix_density * v_int_pt * v_int_pt +
            w * dNdx.transpose() * gamma -
            w * N.transpose() * f * mix_density * std::abs(v_int_pt) *
                v_int_pt / (4 * r_i) -
            w * N.transpose() * Q_mom;
 
        Bh.noalias() +=
            -1 / 2 * w * N.transpose() * rho_dot * v_int_pt * v_int_pt -
            w * N.transpose() * mix_density * v_int_pt * vdot_int_pt +
            1 / 2 * w * dNdx.transpose() * mix_density * v_int_pt * v_int_pt *
                v_int_pt +
            w * N.transpose() * (Q_hx / pi / r_i / r_i) -
            w * N.transpose() * Q_ene;
 
        if (_process_data.has_gravity)
        {
            NodalVectorType gravity_operator =
                N.transpose() * b * w * _element_direction[2];
 
            Bv.noalias() += gravity_operator * mix_density;
            Bh.noalias() += gravity_operator * mix_density * v_int_pt;
        }
    }
 
    //     debugging
    //    std::string sep = "\n----------------------------------------\n";
    //    Eigen::IOFormat CleanFmt(6, 0, ", ", "\n", "[", "]");
    //    std::cout << local_M.format(CleanFmt) << sep;
    //    std::cout << local_K.format(CleanFmt) << sep;
    //    std::cout << local_b.format(CleanFmt) << sep;
}

References _element, _element_direction, _integration_method, _ip_data, _process_data, MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::IAPWSIF97Region1::checkStateInRange(), MaterialPropertyLib::computeVapourVoidFraction(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::driftFluxState(), MaterialPropertyLib::enthalpy, MaterialPropertyLib::VariableArray::enthalpy, enthalpy_index, MaterialPropertyLib::Gas, NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::mixtureSlipParameter(), ProcessLib::WellboreSimulator::NUM_NODAL_DOF, PI, pressure_index, MaterialLib::PhysicalConstant::CriticalPoint::PressureWater, MaterialPropertyLib::saturation_density, MaterialPropertyLib::saturation_enthalpy, ParameterLib::SpatialPosition::setElementID(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::steamDryness(), MaterialPropertyLib::temperature, MaterialPropertyLib::VariableArray::temperature, velocity_index, MaterialPropertyLib::viscosity, and MaterialPropertyLib::voidFractionClosureDiagnostics().

computeSecondaryVariableConcrete()

template<typename ShapeFunction, int GlobalDim>

void ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::computeSecondaryVariableConcrete ( double const , double const , Eigen::VectorXd const & , Eigen::VectorXd const &  )

inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 138 of file WellboreSimulatorFEM.h.

    {
        auto const n_integration_points =
            _integration_method.getNumberOfPoints();
        auto const ele_id = _element.getID();
 
        (*_process_data.mesh_prop_density)[ele_id] =
            std::accumulate(_ip_data.begin(), _ip_data.end(), 0.,
                            [](double const s, auto const& ip)
                            { return s + ip.mix_density; }) /
            n_integration_points;
    }

References _element, _integration_method, _ip_data, and _process_data.

getIntPtDryness()

template<typename ShapeFunction, int GlobalDim>

virtual std::vector< double > const & ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::getIntPtDryness ( const double , std::vector< GlobalVector * > const & , std::vector< NumLib::LocalToGlobalIndexMap const * > const & , std::vector< double > & cache ) const

inlineoverrideprotectedvirtual

Implements ProcessLib::WellboreSimulator::WellboreSimulatorLocalAssemblerInterface.

Definition at line 200 of file WellboreSimulatorFEM.h.

    {
        return ProcessLib::getIntegrationPointScalarData(
            _ip_data, &IpData::dryness, cache);
    }

References _ip_data, ProcessLib::WellboreSimulator::IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType >::dryness, and ProcessLib::getIntegrationPointScalarData().

getIntPtLiquidMassFlowRate()

template<typename ShapeFunction, int GlobalDim>

virtual std::vector< double > const & ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::getIntPtLiquidMassFlowRate ( const double , std::vector< GlobalVector * > const & , std::vector< NumLib::LocalToGlobalIndexMap const * > const & , std::vector< double > & cache ) const

inlineoverrideprotectedvirtual

Implements ProcessLib::WellboreSimulator::WellboreSimulatorLocalAssemblerInterface.

Definition at line 180 of file WellboreSimulatorFEM.h.

    {
        return ProcessLib::getIntegrationPointScalarData(
            _ip_data, &IpData::liquid_mass_flow_rate, cache);
    }

References _ip_data, ProcessLib::getIntegrationPointScalarData(), and ProcessLib::WellboreSimulator::IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType >::liquid_mass_flow_rate.

getIntPtTemperature()

template<typename ShapeFunction, int GlobalDim>

virtual std::vector< double > const & ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::getIntPtTemperature ( const double , std::vector< GlobalVector * > const & , std::vector< NumLib::LocalToGlobalIndexMap const * > const & , std::vector< double > & cache ) const

inlineoverrideprotectedvirtual

Implements ProcessLib::WellboreSimulator::WellboreSimulatorLocalAssemblerInterface.

Definition at line 190 of file WellboreSimulatorFEM.h.

    {
        return ProcessLib::getIntegrationPointScalarData(
            _ip_data, &IpData::temperature, cache);
    }

References _ip_data, ProcessLib::getIntegrationPointScalarData(), and ProcessLib::WellboreSimulator::IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType >::temperature.

getIntPtVaporMassFlowRate()

template<typename ShapeFunction, int GlobalDim>

virtual std::vector< double > const & ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::getIntPtVaporMassFlowRate ( const double , std::vector< GlobalVector * > const & , std::vector< NumLib::LocalToGlobalIndexMap const * > const & , std::vector< double > & cache ) const

inlineoverrideprotectedvirtual

Implements ProcessLib::WellboreSimulator::WellboreSimulatorLocalAssemblerInterface.

Definition at line 170 of file WellboreSimulatorFEM.h.

    {
        return ProcessLib::getIntegrationPointScalarData(
            _ip_data, &IpData::vapor_mass_flow_rate, cache);
    }

References _ip_data, ProcessLib::getIntegrationPointScalarData(), and ProcessLib::WellboreSimulator::IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType >::vapor_mass_flow_rate.

getIntPtVaporVolumeFraction()

template<typename ShapeFunction, int GlobalDim>

virtual std::vector< double > const & ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::getIntPtVaporVolumeFraction ( const double , std::vector< GlobalVector * > const & , std::vector< NumLib::LocalToGlobalIndexMap const * > const & , std::vector< double > & cache ) const

inlineoverrideprotectedvirtual

Implements ProcessLib::WellboreSimulator::WellboreSimulatorLocalAssemblerInterface.

Definition at line 210 of file WellboreSimulatorFEM.h.

    {
        return ProcessLib::getIntegrationPointScalarData(
            _ip_data, &IpData::vapor_volume_fraction, cache);
    }

References _ip_data, ProcessLib::getIntegrationPointScalarData(), and ProcessLib::WellboreSimulator::IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType >::vapor_volume_fraction.

getShapeMatrix()

template<typename ShapeFunction, int GlobalDim>

Eigen::Map< const Eigen::RowVectorXd > ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::getShapeMatrix ( const unsigned integration_point ) const

inlineoverridevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 129 of file WellboreSimulatorFEM.h.

    {
        auto const& N = _ip_data[integration_point].N;
 
        // assumes N is stored contiguously in memory
        return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
    }

References _ip_data.

postTimestepConcrete()

template<typename ShapeFunction, int GlobalDim>

void ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::postTimestepConcrete ( Eigen::VectorXd const & , Eigen::VectorXd const & , double const , double const , int const  )

inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 155 of file WellboreSimulatorFEM.h.

    {
        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();
 
        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            _ip_data[ip].pushBackState();
        }
    }

References _integration_method, and _ip_data.

Member Data Documentation

_element

template<typename ShapeFunction, int GlobalDim>

MeshLib::Element const& ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::_element

protected

Definition at line 220 of file WellboreSimulatorFEM.h.

Referenced by WellboreSimulatorFEM(), assemble(), and computeSecondaryVariableConcrete().

_element_direction

template<typename ShapeFunction, int GlobalDim>

Eigen::Vector3d ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::_element_direction

protected

Definition at line 225 of file WellboreSimulatorFEM.h.

Referenced by WellboreSimulatorFEM(), and assemble().

_integration_method

template<typename ShapeFunction, int GlobalDim>

NumLib::GenericIntegrationMethod const& ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::_integration_method

protected

Definition at line 221 of file WellboreSimulatorFEM.h.

Referenced by WellboreSimulatorFEM(), assemble(), computeSecondaryVariableConcrete(), and postTimestepConcrete().

_ip_data

template<typename ShapeFunction, int GlobalDim>

std::vector<IpData, Eigen::aligned_allocator<IpData> > ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::_ip_data

protected

Definition at line 229 of file WellboreSimulatorFEM.h.

Referenced by WellboreSimulatorFEM(), assemble(), computeSecondaryVariableConcrete(), getIntPtDryness(), getIntPtLiquidMassFlowRate(), getIntPtTemperature(), getIntPtVaporMassFlowRate(), getIntPtVaporVolumeFraction(), getShapeMatrix(), and postTimestepConcrete().

_is_axially_symmetric

template<typename ShapeFunction, int GlobalDim>

bool const ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::_is_axially_symmetric

protected

Definition at line 222 of file WellboreSimulatorFEM.h.

Referenced by WellboreSimulatorFEM().

_process_data

template<typename ShapeFunction, int GlobalDim>

WellboreSimulatorProcessData const& ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::_process_data

protected

Definition at line 223 of file WellboreSimulatorFEM.h.

Referenced by WellboreSimulatorFEM(), assemble(), and computeSecondaryVariableConcrete().

enthalpy_index

template<typename ShapeFunction, int GlobalDim>

const int ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::enthalpy_index = 2 * ShapeFunction::NPOINTS

staticprotected

Definition at line 234 of file WellboreSimulatorFEM.h.

Referenced by assemble().

enthalpy_size

template<typename ShapeFunction, int GlobalDim>

const int ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::enthalpy_size = ShapeFunction::NPOINTS

staticprotected

Definition at line 238 of file WellboreSimulatorFEM.h.

pressure_index

template<typename ShapeFunction, int GlobalDim>

const int ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::pressure_index = 0

staticprotected

Definition at line 232 of file WellboreSimulatorFEM.h.

Referenced by assemble().

pressure_size

template<typename ShapeFunction, int GlobalDim>

const int ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::pressure_size = ShapeFunction::NPOINTS

staticprotected

Definition at line 236 of file WellboreSimulatorFEM.h.

velocity_index

template<typename ShapeFunction, int GlobalDim>

const int ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::velocity_index = ShapeFunction::NPOINTS

staticprotected

Definition at line 233 of file WellboreSimulatorFEM.h.

Referenced by assemble().

velocity_size

template<typename ShapeFunction, int GlobalDim>

const int ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::velocity_size = ShapeFunction::NPOINTS

staticprotected

Definition at line 237 of file WellboreSimulatorFEM.h.

---

The documentation for this class was generated from the following files:

• WellboreSimulatorFEM.h
• WellboreSimulatorFEM-impl.h