Detailed Description

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

Definition at line 31 of file WellboreSimulatorFEM.h.

#include <WellboreSimulatorFEM.h>

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

[legend]

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

[legend]

Public Types
using	ResidualVector = Eigen::Matrix<double, jacobian_residual_size, 1>
using	JacobianMatrix
using	UnknownVector = Eigen::Matrix<double, jacobian_residual_size, 1>

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
void	calculateResidual (double const alpha, double const vapor_water_density, double const liquid_water_density, double const v_mix, double const dryness, double const C_0, double const u_gu, ResidualVector &res)
void	calculateJacobian (double const alpha, double const vapor_water_density, double const liquid_water_density, double const v_mix, double const dryness, double const C_0, double const u_gu, JacobianMatrix &Jac)
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

Static Public Attributes
static int const	jacobian_residual_size = 1

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 49 of file WellboreSimulatorFEM.h.

◆ GlobalDimNodalMatrixType

template<typename ShapeFunction, int GlobalDim>

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

private

Initial value:

typename ShapeMatricesType::GlobalDimNodalMatrixType

EigenFixedShapeMatrixPolicy::GlobalDimNodalMatrixType

MatrixType< GlobalDim, ShapeFunction::NPOINTS > GlobalDimNodalMatrixType

Definition ShapeMatrixPolicy.h:75

Definition at line 47 of file WellboreSimulatorFEM.h.

◆ GlobalDimVectorType

template<typename ShapeFunction, int GlobalDim>

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

private

Definition at line 46 of file WellboreSimulatorFEM.h.

◆ IpData

template<typename ShapeFunction, int GlobalDim>

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

protected

Initial value:

IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>

ProcessLib::WellboreSimulator::IntegrationPointData

Definition WellboreSimulatorLocalAssemblerInterface.h:18

Definition at line 268 of file WellboreSimulatorFEM.h.

◆ JacobianMatrix

template<typename ShapeFunction, int GlobalDim>

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

Initial value:

 
        Eigen::Matrix<double, jacobian_residual_size, jacobian_residual_size,
                      Eigen::RowMajor>

Definition at line 132 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 36 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 39 of file WellboreSimulatorFEM.h.

◆ NodalRowVectorType

template<typename ShapeFunction, int GlobalDim>

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

private

Definition at line 44 of file WellboreSimulatorFEM.h.

◆ NodalVectorType

template<typename ShapeFunction, int GlobalDim>

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

private

Definition at line 43 of file WellboreSimulatorFEM.h.

◆ ResidualVector

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::ResidualVector = Eigen::Matrix<double, jacobian_residual_size, 1>

Definition at line 131 of file WellboreSimulatorFEM.h.

◆ ShapeMatrices

template<typename ShapeFunction, int GlobalDim>

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

private

Definition at line 34 of file WellboreSimulatorFEM.h.

◆ ShapeMatricesType

template<typename ShapeFunction, int GlobalDim>

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

private

Definition at line 33 of file WellboreSimulatorFEM.h.

◆ UnknownVector

template<typename ShapeFunction, int GlobalDim>

using ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::UnknownVector = Eigen::Matrix<double, jacobian_residual_size, 1>

Definition at line 135 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 52 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 22 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;
 
        double liquid_water_density =
            liquid_phase
                .property(MaterialPropertyLib::PropertyType::saturation_density)
                .template value<double>(vars, pos, t, dt);
        double const 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);
 
        // TODO add a function to calculate dryness with constrain of 0
        // to 1.
        double const dryness = std::max(
            0., (h_int_pt - h_sat_liq_w) / (h_sat_vap_w - h_sat_liq_w));
        steam_mass_frac = dryness;
 
        double const 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);
        temperature = T_int_pt;
        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.
 
        // profile parameter of drift flux
        double C_0 = 1 + 0.12 * (1 - dryness);
 
        // For the surface tension calculation, see
        // Cooper, J. R., and R. B. Dooley. "IAPWS release on surface
        // tension of ordinary water substance." International Association
        // for the Properties of Water and Steam (1994).
        double const sigma_gl = 0.2358 *
                                std::pow((1 - T_int_pt / 647.096), 1.256) *
                                (1 - 0.625 * (1 - T_int_pt / 647.096));
        // drift flux velocity
        double const u_gu =
            1.18 * (1 - dryness) *
            std::pow((9.81) * sigma_gl *
                         (liquid_water_density - vapour_water_density),
                     0.25) /
            std::pow(liquid_water_density, 0.5);
 
        // solving void fraction of vapor: Rouhani-Axelsson
        double alpha = 0;
        if (dryness != 0)
        {
            // Local Newton solver
            using LocalJacobianMatrix =
                Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
            using LocalResidualVector = Eigen::Matrix<double, 1, 1>;
            using LocalUnknownVector = Eigen::Matrix<double, 1, 1>;
            LocalJacobianMatrix J_loc;
 
            Eigen::PartialPivLU<LocalJacobianMatrix> linear_solver(1);
 
            auto const update_residual = [&](LocalResidualVector& residual)
            {
                calculateResidual(alpha, vapour_water_density,
                                  liquid_water_density, v_int_pt, dryness, C_0,
                                  u_gu, residual);
            };
 
            auto const update_jacobian = [&](LocalJacobianMatrix& jacobian)
            {
                calculateJacobian(
                    alpha, vapour_water_density, liquid_water_density, v_int_pt,
                    dryness, C_0, u_gu,
                    jacobian);  // for solution dependent Jacobians
            };
 
            auto const update_solution =
                [&](LocalUnknownVector const& increment)
            {
                // increment solution vectors
                alpha += increment[0];
            };
 
            const int maximum_iterations(20);
            const double residuum_tolerance(1.e-10);
            const double increment_tolerance(0);
 
            auto newton_solver = NumLib::NewtonRaphson(
                linear_solver, update_jacobian, update_residual,
                update_solution,
                {maximum_iterations, residuum_tolerance, increment_tolerance});
 
            auto const success_iterations = newton_solver.solve(J_loc);
 
            if (!success_iterations)
            {
                WARN(
                    "Attention! Steam void fraction has not been correctly "
                    "calculated!");
            }
        }
 
        vapor_volume_frac = alpha;
 
        if (alpha == 0)
        {
            liquid_water_density =
                liquid_phase
                    .property(MaterialPropertyLib::PropertyType::density)
                    .template value<double>(vars, pos, t, dt);
        }
 
        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
                         : (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);
 
        // Slip parameter between two phases,
        // see Akbar, Somaieh, N. Fathianpour, and Rafid Al Khoury. "A finite
        // element model for high enthalpy two-phase flow in geothermal
        // wellbores." Renewable Energy 94 (2016): 223-236.
        double const gamma =
            alpha * liquid_water_density * vapour_water_density * mix_density /
            (1 - alpha) /
            std::pow((alpha * C_0 * vapour_water_density +
                      (1 - alpha * C_0) * liquid_water_density),
                     2) *
            std::pow((C_0 - 1) * v_int_pt + u_gu, 2);
 
        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 = std::pow((pi * t_d), -0.5) + 0.5 -
                       0.25 * std::pow((t_d / pi), 0.5) + 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;
}

◆ calculateJacobian()

template<typename ShapeFunction, int GlobalDim>

void ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::calculateJacobian	(	double const	alpha,
		double const	vapor_water_density,
		double const	liquid_water_density,
		double const	v_mix,
		double const	dryness,
		double const	C_0,
		double const	u_gu,
		JacobianMatrix &	Jac )

inline

Definition at line 154 of file WellboreSimulatorFEM.h.

    {
        Jac(0) = dryness * liquid_water_density * v_mix *
                     (vapor_water_density - liquid_water_density) -
                 (C_0 * dryness * liquid_water_density +
                  C_0 * (1 - dryness) * vapor_water_density) *
                     (2 * alpha * vapor_water_density +
                      (1 - 2 * alpha) * liquid_water_density) *
                     v_mix -
                 vapor_water_density * liquid_water_density * u_gu;
    }

Referenced by assemble().

◆ calculateResidual()

template<typename ShapeFunction, int GlobalDim>

void ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::calculateResidual	(	double const	alpha,
		double const	vapor_water_density,
		double const	liquid_water_density,
		double const	v_mix,
		double const	dryness,
		double const	C_0,
		double const	u_gu,
		ResidualVector &	res )

inline

Definition at line 137 of file WellboreSimulatorFEM.h.

    {
        double const rho_mix =
            alpha * vapor_water_density + (1 - alpha) * liquid_water_density;
 
        res(0) =
            dryness * liquid_water_density * rho_mix * v_mix -
            alpha * C_0 * dryness * liquid_water_density * rho_mix * v_mix -
            alpha * C_0 * (1 - dryness) * vapor_water_density * rho_mix *
                v_mix -
            alpha * vapor_water_density * liquid_water_density * u_gu;
    }

Referenced by assemble().

◆ 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 179 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 241 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 221 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 231 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 211 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 251 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 170 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 196 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 261 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 266 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 262 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 270 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 263 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 264 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 275 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 279 of file WellboreSimulatorFEM.h.

◆ jacobian_residual_size

template<typename ShapeFunction, int GlobalDim>

int const ProcessLib::WellboreSimulator::WellboreSimulatorFEM< ShapeFunction, GlobalDim >::jacobian_residual_size = 1

static

Definition at line 130 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 273 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 277 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 274 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 278 of file WellboreSimulatorFEM.h.

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

Detailed Description

Public Types

Public Member Functions

Static Public Attributes

Protected Types

Protected Member Functions

Protected Attributes

Static Protected Attributes

Private Types

Member Typedef Documentation

◆ GlobalDimMatrixType

◆ GlobalDimNodalMatrixType

◆ GlobalDimVectorType

◆ IpData

◆ JacobianMatrix

◆ LocalMatrixType

◆ LocalVectorType

◆ NodalRowVectorType

◆ NodalVectorType

◆ ResidualVector

◆ ShapeMatrices

◆ ShapeMatricesType

◆ UnknownVector

Constructor & Destructor Documentation

◆ WellboreSimulatorFEM()

Member Function Documentation

◆ assemble()

◆ calculateJacobian()

◆ calculateResidual()

◆ computeSecondaryVariableConcrete()

◆ getIntPtDryness()

◆ getIntPtLiquidMassFlowRate()

◆ getIntPtTemperature()

◆ getIntPtVaporMassFlowRate()

◆ getIntPtVaporVolumeFraction()

◆ getShapeMatrix()

◆ postTimestepConcrete()

Member Data Documentation

◆ _element

◆ _element_direction

◆ _integration_method

◆ _ip_data

◆ _is_axially_symmetric

◆ _process_data

◆ enthalpy_index

◆ enthalpy_size

◆ jacobian_residual_size

◆ pressure_index

◆ pressure_size

◆ velocity_index

◆ velocity_size