OGS
ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > Class Template Reference

Detailed Description

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
class ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >

Definition at line 45 of file ThermoHydroMechanicsFEM.h.

#include <ThermoHydroMechanicsFEM.h>

Inheritance diagram for ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:
[legend]
Collaboration diagram for ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:
[legend]

Public Types

using ShapeMatricesTypeDisplacement
using ShapeMatricesTypePressure
using GlobalDimMatrixType
using GlobalDimVectorType
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>
using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>

Public Member Functions

 ThermoHydroMechanicsLocalAssembler (ThermoHydroMechanicsLocalAssembler const &)=delete
 ThermoHydroMechanicsLocalAssembler (ThermoHydroMechanicsLocalAssembler &&)=delete
 ThermoHydroMechanicsLocalAssembler (MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, ThermoHydroMechanicsProcessData< DisplacementDim > &process_data)
std::size_t setIPDataInitialConditions (std::string_view const name, double const *values, int const integration_order) override
 Returns number of read integration points.
void assemble (double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
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_rhs_data, std::vector< double > &local_Jac_data) override
void initializeConcrete () override
void setInitialConditionsConcrete (Eigen::VectorXd const local_x, double const t, int const process_id) override
void preTimestepConcrete (std::vector< double > const &, double const, double const) override
void postTimestepConcrete (Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, int const) override
void computeSecondaryVariableConcrete (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev) override
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const override
 Provides the shape matrix at the given integration point.
std::vector< double > const & getIntPtDarcyVelocity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double > getSigma () const override
std::vector< double > getSigmaIce () const override
std::vector< double > getIceVolumeFraction () const override
std::vector< double > const & getIntPtFluidDensity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double > const & getIntPtViscosity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
int getNumberOfVectorElementsForDeformation () const override
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 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.
Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual ~ExtrapolatableElement ()=default

Static Public Attributes

static int const KelvinVectorSize
static constexpr auto & N_u_op

Private Types

using BMatricesType
using IpData

Private Member Functions

ConstitutiveRelationsValues< DisplacementDim > updateConstitutiveRelations (Eigen::Ref< Eigen::VectorXd const > const local_x, Eigen::Ref< Eigen::VectorXd const > const local_x_prev, ParameterLib::SpatialPosition const &x_position, double const t, double const dt, IpData &ip_data, IntegrationPointDataForOutput< DisplacementDim > &ip_data_output) const
std::size_t setSigma (double const *values)
std::vector< double > const & getIntPtSigma (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
std::vector< double > const & getIntPtSigmaIce (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
std::vector< double > getEpsilon0 () const override
virtual std::vector< double > const & getIntPtEpsilon0 (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
std::vector< double > getEpsilonM () const override
virtual std::vector< double > const & getIntPtEpsilonM (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
std::vector< double > getEpsilon () const override
virtual std::vector< double > const & getIntPtEpsilon (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
std::vector< double > const & getIntPtIceVolume (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
unsigned getNumberOfIntegrationPoints () const override
int getMaterialID () const override
std::vector< double > getMaterialStateVariableInternalState (std::function< std::span< double >(typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables &)> const &get_values_span, int const &n_components) const override
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt (unsigned integration_point) const override

Static Private Member Functions

template<typename SolutionVector>
static constexpr auto localDOF (SolutionVector const &x)

Private Attributes

ThermoHydroMechanicsProcessData< DisplacementDim > & _process_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
std::vector< IntegrationPointDataForOutput< DisplacementDim >, Eigen::aligned_allocator< IntegrationPointDataForOutput< DisplacementDim > > > _ip_data_output
NumLib::GenericIntegrationMethod const & _integration_method
MeshLib::Element const & _element
bool const _is_axially_symmetric
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data

Static Private Attributes

static const int temperature_index = 0
static const int temperature_size = ShapeFunctionPressure::NPOINTS
static const int pressure_index = ShapeFunctionPressure::NPOINTS
static const int pressure_size = ShapeFunctionPressure::NPOINTS
static const int displacement_index = ShapeFunctionPressure::NPOINTS * 2
static const int displacement_size

Member Typedef Documentation

◆ BMatricesType

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
using ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::BMatricesType
private

◆ GlobalDimMatrixType

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
using ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::GlobalDimMatrixType
Initial value:

Definition at line 56 of file ThermoHydroMechanicsFEM.h.

◆ GlobalDimVectorType

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
using ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::GlobalDimVectorType
Initial value:

Definition at line 59 of file ThermoHydroMechanicsFEM.h.

◆ Invariants

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
using ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>

Definition at line 64 of file ThermoHydroMechanicsFEM.h.

◆ IpData

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
using ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::IpData
private
Initial value:
ShapeMatricesTypePressure, DisplacementDim,
ShapeFunctionDisplacement::NPOINTS>
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > BMatricesType

Definition at line 286 of file ThermoHydroMechanicsFEM.h.

◆ ShapeMatricesTypeDisplacement

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
using ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypeDisplacement
Initial value:
EigenFixedShapeMatrixPolicy< ShapeFunction, GlobalDim > ShapeMatrixPolicyType

Definition at line 49 of file ThermoHydroMechanicsFEM.h.

◆ ShapeMatricesTypePressure

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
using ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypePressure

◆ SymmetricTensor

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
using ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>

Definition at line 66 of file ThermoHydroMechanicsFEM.h.

Constructor & Destructor Documentation

◆ ThermoHydroMechanicsLocalAssembler() [1/3]

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ThermoHydroMechanicsLocalAssembler ( ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > const & )
delete

◆ ThermoHydroMechanicsLocalAssembler() [2/3]

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ThermoHydroMechanicsLocalAssembler ( ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > && )
delete

◆ ThermoHydroMechanicsLocalAssembler() [3/3]

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ThermoHydroMechanicsLocalAssembler ( MeshLib::Element const & e,
std::size_t const ,
NumLib::GenericIntegrationMethod const & integration_method,
bool const is_axially_symmetric,
ThermoHydroMechanicsProcessData< DisplacementDim > & process_data )

Definition at line 32 of file ThermoHydroMechanicsFEM-impl.h.

41 _element(e),
43{
44 unsigned const n_integration_points =
45 _integration_method.getNumberOfPoints();
46
50
51 auto const shape_matrices_u =
56
57 auto const shape_matrices_p =
61
62 auto const& solid_material =
64 _process_data.solid_materials, _process_data.material_ids,
65 e.getID());
66
67 // Consistency check: if frozen liquid phase is given, then the constitutive
68 // relation for ice must also be given, and vice versa.
69 auto const& medium = _process_data.media_map.getMedium(_element.getID());
71 (_process_data.ice_constitutive_relation != nullptr))
72 {
74 "Frozen liquid phase is {:s} and the solid material constitutive "
75 "relation for ice is {:s}. But both must be given (or both "
76 "omitted).",
78 ? "specified"
79 : "not specified",
80 _process_data.ice_constitutive_relation != nullptr
81 ? "specified"
82 : "not specified");
83 }
84 for (unsigned ip = 0; ip < n_integration_points; ip++)
85 {
86 _ip_data.emplace_back(solid_material);
87 auto& ip_data = _ip_data[ip];
88 auto const& sm_u = shape_matrices_u[ip];
89 ip_data.integration_weight =
90 _integration_method.getWeightedPoint(ip).getWeight() *
91 sm_u.integralMeasure * sm_u.detJ;
92
93 ip_data.N_u = sm_u.N;
94 ip_data.dNdx_u = sm_u.dNdx;
95
97 ip_data.dNdx = shape_matrices_p[ip].dNdx;
98
100 }
101}
#define OGS_FATAL(...)
Definition Error.h:10
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
ThermoHydroMechanicsProcessData< DisplacementDim > & _process_data
std::vector< IntegrationPointDataForOutput< DisplacementDim >, Eigen::aligned_allocator< IntegrationPointDataForOutput< DisplacementDim > > > _ip_data_output
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)

References _element, _integration_method, _ip_data, _ip_data_output, _is_axially_symmetric, _process_data, _secondary_data, MaterialPropertyLib::FrozenLiquid, MeshLib::Element::getID(), NumLib::initShapeMatrices(), OGS_FATAL, and MaterialLib::Solids::selectSolidConstitutiveRelation().

Member Function Documentation

◆ assemble()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
void ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assemble ( double const ,
double const ,
std::vector< double > const & ,
std::vector< double > const & ,
std::vector< double > & ,
std::vector< double > & ,
std::vector< double > &  )
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 90 of file ThermoHydroMechanicsFEM.h.

96 {
98 "ThermoHydroMechanicsLocalAssembler: assembly without Jacobian is "
99 "not implemented.");
100 }

References OGS_FATAL.

◆ assembleWithJacobian()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
void ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobian ( double const t,
double const dt,
std::vector< double > const & local_x,
std::vector< double > const & local_x_prev,
std::vector< double > & local_rhs_data,
std::vector< double > & local_Jac_data )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 696 of file ThermoHydroMechanicsFEM-impl.h.

702{
703 assert(local_x.size() ==
705
706 auto const x =
709 local_x_prev.size());
710
711 auto const [T, p, u] = localDOF(local_x);
712 auto const [T_prev, p_prev, u_prev] = localDOF(local_x_prev);
713
720
725
728
731
734
737
740
743
746
749
752
755 Kup;
757
760 Kpu;
761 if (!_process_data.is_volume_balance_equation_type)
762 {
764 }
765
766 auto const& medium = _process_data.media_map.getMedium(_element.getID());
767 bool const has_frozen_liquid_phase =
769
770 unsigned const n_integration_points =
771 _integration_method.getNumberOfPoints();
772
774 double average_velocity_norm = 0.0;
776
777 for (unsigned ip = 0; ip < n_integration_points; ip++)
778 {
779 auto& ip_data = _ip_data[ip];
780 auto const& N_u = ip_data.N_u;
782 std::nullopt, _element.getID(),
786 _element, N_u))};
787
788 auto const crv = updateConstitutiveRelations(
790
791 auto const& w = ip_data.integration_weight;
792
793 auto const& dNdx_u = ip_data.dNdx_u;
794
795 auto const& N = ip_data.N;
796 auto const& dNdx = ip_data.dNdx;
797
798 auto const T_int_pt = N.dot(T);
799
800 auto const x_coord =
801 x_position.getCoordinates().value()[0]; // r for axisymmetry
802 auto const B =
807
808 auto const& b = _process_data.specific_body_force;
809 auto const velocity = _ip_data_output[ip].velocity;
810
811 //
812 // displacement equation, displacement part
813 //
814
815 auto const C_eff = has_frozen_liquid_phase
816 ? (crv.C + crv.J_uu_fr).eval()
817 : crv.C.eval();
821 .noalias() += B.transpose() * C_eff * B * w;
822
823 auto const uT_coeff =
825 ? (crv.J_uT_fr +
826 crv.C * crv.solid_linear_thermal_expansion_coefficient)
827 .eval()
829 .eval();
830
832 {
834 .noalias() -= B.transpose() * crv.r_u_fr * w;
835 }
836
840 .noalias() -= B.transpose() * uT_coeff * N * w;
841
843 .noalias() -= (B.transpose() * ip_data.sigma_eff -
844 N_u_op(N_u).transpose() * crv.rho * b) *
845 w;
846
847 //
848 // displacement equation, pressure part (K_up)
849 //
850 double const fluid_density = _ip_data_output[ip].fluid_density;
851 double const up_coeff =
853 ? crv.alpha_biot * ip_data.phi_fr / ip_data.porosity *
854 (_ip_data_output[ip].rho_fr / fluid_density - 1) +
855 crv.alpha_biot
856 : crv.alpha_biot;
857 Kup.noalias() +=
858 B.transpose() * Invariants::identity2 * N * (up_coeff * w);
859
860 //
861 // pressure equation, pressure part (K_pp and M_pp).
862 //
863 double const scaling_factor =
864 _process_data.is_volume_balance_equation_type ? 1.0 : fluid_density;
865
866 laplace_p.noalias() += dNdx.transpose() * crv.K_over_mu * dNdx *
867 (crv.k_rel * scaling_factor * w);
871 .noalias() += dNdx.transpose() * crv.K_over_mu * (dNdx * p) * N *
872 (crv.dk_rel_dT * scaling_factor * w);
873
874 double const storage_p_coeff_no_fr =
875 ip_data.porosity * crv.fluid_compressibility +
876 (crv.alpha_biot - ip_data.porosity) * crv.beta_SR;
877 double const storage_p_coeff =
880
881 storage_p.noalias() +=
882 N.transpose() * N * (storage_p_coeff * scaling_factor * w);
883
885 {
889 .noalias() +=
890 N.transpose() * crv.J_pT_fr * N * scaling_factor * w;
891 }
892
893 laplace_T.noalias() += dNdx.transpose() * crv.K_pT_thermal_osmosis *
895 //
896 // RHS, pressure part
897 //
898 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
899 N * (up_coeff * crv.eps_v_dot * scaling_factor * w);
900
901 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
902 dNdx.transpose() * crv.K_over_mu * b *
903 (fluid_density * crv.k_rel * scaling_factor * w);
907 .noalias() -= dNdx.transpose() * crv.K_over_mu * b * N *
908 (fluid_density * crv.dk_rel_dT * scaling_factor * w);
909
910 //
911 // pressure equation, temperature part (M_pT)
912 //
913
914 double const storage_T_coeff =
915 has_frozen_liquid_phase ? crv.storage_T_fr + crv.beta : crv.beta;
916
917 storage_T.noalias() +=
918 N.transpose() * storage_T_coeff * N * scaling_factor * w;
919
920 //
921 // pressure equation, displacement part.
922 //
923 // reusing Kup.transpose() if the equation balance type is not volume.
924
925 if (!_process_data.is_volume_balance_equation_type)
926 {
927 Kpu.noalias() += N.transpose() * Invariants::identity2.transpose() *
928 B * (up_coeff * scaling_factor * w);
929
930 //
931 // The contribution to Jacobian from d()/ drho
932 // drho/dp, d()/ drho drho/dT:
933 //
934 double const storage_p_solid_coeff =
935 (crv.alpha_biot - ip_data.porosity) * crv.beta_SR;
936
937 double const p_dot = N.dot(p - p_prev) / dt;
938 double const T_dot = N.dot(T - T_prev) / dt;
939 double const drho_dp_coeff = storage_p_solid_coeff * p_dot +
941 up_coeff * crv.eps_v_dot;
942
946 .noalias() +=
947 // TODO (WW) : Add ip_data.porosity * d2rho_LR_dp2 * w.
948 N.transpose() * N * (drho_dp_coeff * crv.drho_LR_dp * w);
952 .noalias() +=
953 // TODO (WW) : Add ip_data.porosity * d2rho_LR_dpdT * w.
954 N.transpose() * N * (drho_dp_coeff * crv.drho_LR_dT * w);
955
956 // The term from d (rho_L K(grad p - rho_L b)/dp:
957 // derivative of rhp_L * K_over_mu * k_rel (grad p + rho_l g) with
958 // respect to pressure and temperature.
959 auto const dlaplace_temporal_factor =
960 (-velocity - (fluid_density * crv.k_rel) * crv.K_over_mu * b);
964 .noalias() += dNdx.transpose() * dlaplace_temporal_factor * N *
965 (crv.drho_LR_dp * w);
969 .noalias() += dNdx.transpose() * dlaplace_temporal_factor * N *
970 (crv.drho_LR_dT * w);
971 }
972
973 //
974 // temperature equation, temperature part.
975 //
976 KTT.noalias() +=
977 dNdx.transpose() * crv.effective_thermal_conductivity * dNdx * w;
978 dKTT_dT_T.noalias() +=
979 dNdx.transpose() * crv.dlambda_eff_dT * dNdx * T * N * w;
980
981 ip_flux_vector.emplace_back(velocity * fluid_density * crv.c_f);
982 // Without any flux correction the flux derivative is as follows. The
983 // contribution to KTT is different if any stabilization scheme is used,
984 // but this is ignored for the moment.
986 crv.dvelocity_dT * fluid_density * crv.c_f +
987 velocity * crv.drho_LR_dT * crv.c_f;
988 dKTT_dT_T.noalias() +=
989 N.transpose() * dip_flux_vector_dT.transpose() * dNdx * T * N * w;
991
992 MTT.noalias() += N.transpose() *
993 (crv.sensible_volumetric_heat_capacity -
994 crv.latent_volumetric_heat_capacity) *
995 N * w;
999 .noalias() += N.transpose() * crv.J_TT * N * w;
1000
1001 //
1002 // temperature equation, pressure part
1003 //
1004 KTp.noalias() +=
1005 dNdx.transpose() * crv.K_pT_thermal_osmosis * dNdx * (T_int_pt * w);
1006
1007 // linearized darcy
1008 dKTT_dp_T.noalias() -= N.transpose() * (dNdx * T).transpose() *
1009 crv.K_over_mu * dNdx *
1010 (fluid_density * crv.c_f * crv.k_rel * w);
1011
1012 /* TODO (Joerg) Temperature changes due to thermal dilatation of the
1013 * fluid, which are usually discarded as being very small.
1014 * Zhou et al. (10.1016/S0020-7683(98)00089-4) states that:
1015 * "Biot (1956) neglected this term and it is included here for
1016 * completeness"
1017 * Keeping the code here in the case these are needed for the named
1018 * effects in the future.
1019 if (fluid_compressibility != 0)
1020 {
1021 auto const C_el = ip_data.computeElasticTangentStiffness(
1022 t, x_position, dt, static_cast<double>(T_int_pt));
1023 auto const solid_skeleton_compressibility =
1024 1 / solid_material.getBulkModulus(t, x_position, &C_el);
1025 double const fluid_volumetric_thermal_expansion_coefficient =
1026 MaterialPropertyLib::getLiquidThermalExpansivity(
1027 liquid_phase, vars, fluid_density, x_position, t, dt);
1028
1029 KTT.noalias() +=
1030 dNdx.transpose() *
1031 (-T_int_pt * fluid_volumetric_thermal_expansion_coefficient *
1032 K_pT_thermal_osmosis / fluid_compressibility) *
1033 dNdx * w;
1034
1035 local_rhs.template segment<temperature_size>(temperature_index)
1036 .noalias() +=
1037 dNdx.transpose() *
1038 (-T_int_pt * fluid_volumetric_thermal_expansion_coefficient /
1039 fluid_compressibility) *
1040 fluid_density * crv.k_rel * K_over_mu * b * w;
1041 MTu part for rhs and Jacobian:
1042 (-T_int_pt *
1043 Invariants::trace(solid_linear_thermal_expansion_coefficient) /
1044 solid_skeleton_compressibility) *
1045 N.transpose() * identity2.transpose() * B * w;
1046 KTp part for rhs and Jacobian:
1047 dNdx.transpose() *
1048 (T_int_pt * fluid_volumetric_thermal_expansion_coefficient *
1049 crv.k_rel * K_over_mu / fluid_compressibility) *
1050 dNdx * w;
1051 }
1052 */
1053 }
1054
1057 average_velocity_norm / static_cast<double>(n_integration_points), KTT);
1058
1059 // temperature equation, temperature part
1060 local_Jac
1063 .noalias() += KTT + dKTT_dT_T + MTT / dt;
1064
1065 // temperature equation, pressure part
1066 local_Jac
1069 .noalias() += KTp + dKTT_dp_T;
1070
1071 // displacement equation, pressure part
1072 local_Jac
1075 .noalias() -= Kup;
1076
1077 // pressure equation, temperature part.
1078 local_Jac
1081 .noalias() += -storage_T / dt + laplace_T;
1082
1083 // pressure equation, pressure part.
1084 local_Jac
1087 .noalias() += laplace_p + storage_p / dt;
1088
1089 // pressure equation, displacement part.
1090 if (_process_data.is_volume_balance_equation_type)
1091 {
1092 local_Jac
1095 .noalias() += Kup.transpose() / dt;
1096 }
1097 else
1098 {
1099 local_Jac
1102 .noalias() += Kpu / dt;
1103 }
1104
1105 // pressure equation (f_p)
1106 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
1107 laplace_p * p + laplace_T * T + storage_p * (p - p_prev) / dt -
1108 storage_T * (T - T_prev) / dt;
1109
1110 // displacement equation (f_u)
1112 .noalias() += Kup * p;
1113
1114 // temperature equation (f_T)
1116 KTT * T + MTT * (T - T_prev) / dt;
1117
1119 KTp * p;
1120}
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
ConstitutiveRelationsValues< DisplacementDim > updateConstitutiveRelations(Eigen::Ref< Eigen::VectorXd const > const local_x, Eigen::Ref< Eigen::VectorXd const > const local_x_prev, ParameterLib::SpatialPosition const &x_position, double const t, double const dt, IpData &ip_data, IntegrationPointDataForOutput< DisplacementDim > &ip_data_output) const
typename ShapeMatricesTypePressure::GlobalDimVectorType GlobalDimVectorType
void assembleAdvectionMatrix(IPData const &ip_data_vector, NumLib::ShapeMatrixCache const &shape_matrix_cache, std::vector< FluxVectorType > const &ip_flux_vector, Eigen::MatrixBase< Derived > &laplacian_matrix)
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.

References _element, _integration_method, _ip_data, _ip_data_output, _is_axially_symmetric, _process_data, NumLib::detail::assembleAdvectionMatrix(), assembleWithJacobian(), ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), displacement_index, displacement_size, MaterialPropertyLib::FrozenLiquid, ParameterLib::SpatialPosition::getCoordinates(), MathLib::KelvinVector::Invariants< KelvinVectorSize >::identity2, NumLib::interpolateCoordinates(), localDOF(), N_u_op, pressure_index, pressure_size, temperature_index, temperature_size, and updateConstitutiveRelations().

Referenced by assembleWithJacobian().

◆ computeSecondaryVariableConcrete()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
void ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::computeSecondaryVariableConcrete ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1181 of file ThermoHydroMechanicsFEM-impl.h.

1185{
1186 auto const p = local_x.template segment<pressure_size>(pressure_index);
1187 auto const T =
1189
1190 unsigned const n_integration_points =
1191 _integration_method.getNumberOfPoints();
1192
1193 double phi_fr_avg = 0;
1194 double fluid_density_avg = 0;
1195 double viscosity_avg = 0;
1196
1198 KV sigma_avg = KV::Zero();
1200
1201 for (unsigned ip = 0; ip < n_integration_points; ip++)
1202 {
1203 auto& ip_data = _ip_data[ip];
1204
1205 phi_fr_avg += ip_data.phi_fr;
1206 fluid_density_avg += _ip_data_output[ip].fluid_density;
1207 viscosity_avg += _ip_data_output[ip].viscosity;
1208 sigma_avg += ip_data.sigma_eff;
1209 sigma_ice_avg += ip_data.sigma_eff_ice;
1210 }
1211
1217
1218 (*_process_data.element_phi_fr)[_element.getID()] = phi_fr_avg;
1219 (*_process_data.element_fluid_density)[_element.getID()] =
1221 (*_process_data.element_viscosity)[_element.getID()] = viscosity_avg;
1222
1223 Eigen::Map<KV>(&(*_process_data.element_stresses)[_element.getID() *
1226
1231
1235 *_process_data.pressure_interpolated);
1236
1240 *_process_data.temperature_interpolated);
1241}
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)

References _element, _integration_method, _ip_data, _ip_data_output, _is_axially_symmetric, _process_data, NumLib::interpolateToHigherOrderNodes(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), pressure_index, and temperature_index.

◆ getEpsilon()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
std::vector< double > ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getEpsilon ( ) const
inlineoverrideprivatevirtual

Implements ProcessLib::ThermoHydroMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 364 of file ThermoHydroMechanicsFEM.h.

365 {
366 constexpr int kelvin_vector_size =
368
371 { return getIntPtEpsilon(0, {}, {}, values); });
372 }
virtual std::vector< double > const & getIntPtEpsilon(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
std::vector< double > transposeInPlace(StoreValuesFunction const &store_values_function)

References getIntPtEpsilon(), MathLib::KelvinVector::kelvin_vector_dimensions(), and ProcessLib::transposeInPlace().

◆ getEpsilon0()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
std::vector< double > ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getEpsilon0 ( ) const
inlineoverrideprivatevirtual

Implements ProcessLib::ThermoHydroMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 325 of file ThermoHydroMechanicsFEM.h.

326 {
327 constexpr int kelvin_vector_size =
329
332 { return getIntPtEpsilon0(0, {}, {}, values); });
333 }
virtual std::vector< double > const & getIntPtEpsilon0(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override

References getIntPtEpsilon0(), MathLib::KelvinVector::kelvin_vector_dimensions(), and ProcessLib::transposeInPlace().

◆ getEpsilonM()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
std::vector< double > ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getEpsilonM ( ) const
inlineoverrideprivatevirtual

Implements ProcessLib::ThermoHydroMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 344 of file ThermoHydroMechanicsFEM.h.

345 {
346 constexpr int kelvin_vector_size =
348
351 { return getIntPtEpsilonM(0, {}, {}, values); });
352 }
virtual std::vector< double > const & getIntPtEpsilonM(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override

References getIntPtEpsilonM(), MathLib::KelvinVector::kelvin_vector_dimensions(), and ProcessLib::transposeInPlace().

◆ getIceVolumeFraction()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
std::vector< double > ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIceVolumeFraction ( ) const
inlineoverridevirtual

Implements ProcessLib::ThermoHydroMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 259 of file ThermoHydroMechanicsFEM.h.

260 {
262 getIntPtIceVolume(0, {}, {}, result);
263 return result;
264 }
std::vector< double > const & getIntPtIceVolume(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override

References getIntPtIceVolume().

◆ getIntPtDarcyVelocity()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
std::vector< double > const & ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtDarcyVelocity ( const double t,
std::vector< GlobalVector * > const & x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table,
std::vector< double > & cache ) const
overridevirtual

Implements ProcessLib::ThermoHydroMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1125 of file ThermoHydroMechanicsFEM-impl.h.

1131{
1132 unsigned const n_integration_points =
1133 _integration_method.getNumberOfPoints();
1134
1135 cache.clear();
1139
1140 for (unsigned ip = 0; ip < n_integration_points; ip++)
1141 {
1142 cache_matrix.col(ip).noalias() = _ip_data_output[ip].velocity;
1143 }
1144
1145 return cache;
1146}

References _integration_method, _ip_data_output, and MathLib::createZeroedMatrix().

◆ getIntPtEpsilon()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
virtual std::vector< double > const & ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEpsilon ( const double ,
std::vector< GlobalVector * > const & ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & ,
std::vector< double > & cache ) const
inlineoverrideprivatevirtual

◆ getIntPtEpsilon0()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
virtual std::vector< double > const & ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEpsilon0 ( const double ,
std::vector< GlobalVector * > const & ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & ,
std::vector< double > & cache ) const
inlineoverrideprivatevirtual

◆ getIntPtEpsilonM()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
virtual std::vector< double > const & ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEpsilonM ( const double ,
std::vector< GlobalVector * > const & ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & ,
std::vector< double > & cache ) const
inlineoverrideprivatevirtual

◆ getIntPtFluidDensity()

template<typename ShapeFunctionDisplacement, typename ShapeFunction, int DisplacementDim>
std::vector< double > const & ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, DisplacementDim >::getIntPtFluidDensity ( const double t,
std::vector< GlobalVector * > const & x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table,
std::vector< double > & cache ) const
overridevirtual

◆ getIntPtIceVolume()

◆ getIntPtSigma()

◆ getIntPtSigmaIce()

◆ getIntPtViscosity()

template<typename ShapeFunctionDisplacement, typename ShapeFunction, int DisplacementDim>
std::vector< double > const & ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunction, DisplacementDim >::getIntPtViscosity ( const double t,
std::vector< GlobalVector * > const & x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table,
std::vector< double > & cache ) const
overridevirtual

◆ getMaterialID()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
int ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMaterialID ( ) const
inlineoverrideprivatevirtual

Implements ProcessLib::ThermoHydroMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 399 of file ThermoHydroMechanicsFEM.h.

400 {
401 return _process_data.material_ids == nullptr
402 ? 0
403 : (*_process_data.material_ids)[_element.getID()];
404 }

References _element, and _process_data.

◆ getMaterialStateVariableInternalState()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
std::vector< double > ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMaterialStateVariableInternalState ( std::function< std::span< double >(typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables &)> const & get_values_span,
int const & n_components ) const
inlineoverrideprivatevirtual

◆ getMaterialStateVariablesAt()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMaterialStateVariablesAt ( unsigned integration_point) const
inlineoverrideprivatevirtual

Implements ProcessLib::ThermoHydroMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 419 of file ThermoHydroMechanicsFEM.h.

420 {
421 return *_ip_data[integration_point].material_state_variables;
422 }

References _ip_data.

◆ getNumberOfIntegrationPoints()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
unsigned ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getNumberOfIntegrationPoints ( ) const
inlineoverrideprivatevirtual

Implements ProcessLib::ThermoHydroMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 394 of file ThermoHydroMechanicsFEM.h.

395 {
396 return _integration_method.getNumberOfPoints();
397 }

References _integration_method.

◆ getNumberOfVectorElementsForDeformation()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
int ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getNumberOfVectorElementsForDeformation ( ) const
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 278 of file ThermoHydroMechanicsFEM.h.

279 {
280 return displacement_size;
281 }

References displacement_size.

◆ getShapeMatrix()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
Eigen::Map< const Eigen::RowVectorXd > ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getShapeMatrix ( const unsigned integration_point) const
inlineoverridevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 221 of file ThermoHydroMechanicsFEM.h.

223 {
224 auto const& N_u = _secondary_data.N_u[integration_point];
225
226 // assumes N is stored contiguously in memory
227 return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
228 }

References _secondary_data.

◆ getSigma()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
std::vector< double > ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getSigma ( ) const
inlineoverridevirtual

Implements ProcessLib::ThermoHydroMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 239 of file ThermoHydroMechanicsFEM.h.

240 {
241 constexpr int kelvin_vector_size =
243
246 { return getIntPtSigma(0, {}, {}, values); });
247 }
std::vector< double > const & getIntPtSigma(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override

References getIntPtSigma(), MathLib::KelvinVector::kelvin_vector_dimensions(), and ProcessLib::transposeInPlace().

◆ getSigmaIce()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
std::vector< double > ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getSigmaIce ( ) const
inlineoverridevirtual

Implements ProcessLib::ThermoHydroMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 249 of file ThermoHydroMechanicsFEM.h.

250 {
251 constexpr int kelvin_vector_size =
253
256 { return getIntPtSigmaIce(0, {}, {}, values); });
257 }
std::vector< double > const & getIntPtSigmaIce(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override

References getIntPtSigmaIce(), MathLib::KelvinVector::kelvin_vector_dimensions(), and ProcessLib::transposeInPlace().

◆ initializeConcrete()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
void ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete ( )
inlineoverridevirtual

Set initial stress from parameter.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 108 of file ThermoHydroMechanicsFEM.h.

109 {
110 unsigned const n_integration_points =
111 _integration_method.getNumberOfPoints();
112
113 for (unsigned ip = 0; ip < n_integration_points; ip++)
114 {
115 auto& ip_data = _ip_data[ip];
116
118 std::nullopt, _element.getID(),
123
125 if (_process_data.initial_stress.value)
126 {
127 ip_data.sigma_eff =
129 DisplacementDim>((*_process_data.initial_stress.value)(
131 independent
132 */
133 ,
134 x_position));
135 }
136
137 double const t = 0; // TODO (naumov) pass t from top
138 ip_data.solid_material.initializeInternalStateVariables(
139 t, x_position, *ip_data.material_state_variables);
140
141 ip_data.pushBackState();
142 }
143 }

References _element, _integration_method, _ip_data, _process_data, NumLib::interpolateCoordinates(), and MathLib::KelvinVector::symmetricTensorToKelvinVector().

◆ localDOF()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
template<typename SolutionVector>
constexpr auto ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::localDOF ( SolutionVector const & x)
inlinestaticconstexprprivate

◆ postTimestepConcrete()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
void ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::postTimestepConcrete ( Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev,
double const t,
double const dt,
int const  )
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 162 of file ThermoHydroMechanicsFEM.h.

166 {
167 unsigned const n_integration_points =
168 _integration_method.getNumberOfPoints();
169
170 auto const [T, p, u] = localDOF(local_x);
171 auto const [T_prev, p_prev, u_prev] = localDOF(local_x_prev);
172
173 // Stored sensible heat integrated over the element, for output.
174 double sensible_heat = 0.0;
175
176 for (unsigned ip = 0; ip < n_integration_points; ip++)
177 {
178 auto& ip_data = _ip_data[ip];
179 auto const& N_u = ip_data.N_u;
180 auto const& dNdx_u = ip_data.dNdx_u;
181
183 std::nullopt, _element.getID(),
188
189 auto const crv = updateConstitutiveRelations(
192
193 // Stored sensible heat (latent already excluded).
194 sensible_heat += crv.sensible_volumetric_heat_capacity *
195 ip_data.N.dot(T) * ip_data.integration_weight;
196
197 auto const x_coord =
198 x_position.getCoordinates().value()[0]; // r for axisymmetry
199 auto const B = LinearBMatrix::computeBMatrix<
203
205 eps_prev = B * u_prev;
206
207 _ip_data[ip].eps0 =
208 _ip_data[ip].eps0_prev +
209 (1 - _ip_data[ip].phi_fr_prev / _ip_data[ip].porosity) *
211 _ip_data[ip].pushBackState();
212 }
213
214 (*_process_data.cell_sensible_heat)[_element.getID()] = sensible_heat;
215 }

References _element, _integration_method, _ip_data, _ip_data_output, _is_axially_symmetric, _process_data, ProcessLib::LinearBMatrix::computeBMatrix(), ParameterLib::SpatialPosition::getCoordinates(), NumLib::interpolateCoordinates(), localDOF(), and updateConstitutiveRelations().

◆ preTimestepConcrete()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
void ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::preTimestepConcrete ( std::vector< double > const & ,
double const ,
double const  )
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 149 of file ThermoHydroMechanicsFEM.h.

151 {
152 unsigned const n_integration_points =
153 _integration_method.getNumberOfPoints();
154
155 for (unsigned ip = 0; ip < n_integration_points; ip++)
156 {
157 _ip_data_output[ip].velocity.setConstant(
159 }
160 }

References _integration_method, and _ip_data_output.

◆ setInitialConditionsConcrete()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
void ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::setInitialConditionsConcrete ( Eigen::VectorXd const local_x,
double const t,
int const process_id )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 177 of file ThermoHydroMechanicsFEM-impl.h.

181{
182 // TODO: For staggered scheme, overload
183 // LocalAssemblerInterface::setInitialConditions to enable local_x contains
184 // the primary variables from all coupled processes.
185 auto const [T, p, u] = localDOF(local_x);
186
187 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
188
190 auto const& medium = _process_data.media_map.getMedium(_element.getID());
191 auto* const frozen_liquid_phase =
194 : nullptr;
195
196 int const n_integration_points = _integration_method.getNumberOfPoints();
197 for (int ip = 0; ip < n_integration_points; ip++)
198 {
199 auto& ip_data = _ip_data[ip];
200 auto const& N = ip_data.N;
201 auto const& N_u = ip_data.N_u;
203 std::nullopt, _element.getID(),
207 _element, N_u))};
208
209 auto& sigma_eff = ip_data.sigma_eff;
210 if (_process_data.initial_stress.isTotalStress())
211 {
212 auto const alpha_b =
214 .template value<double>(vars, x_position, t, dt);
215
216 sigma_eff.noalias() += alpha_b * N.dot(p) * Invariants::identity2;
217 }
218 ip_data.sigma_eff_prev.noalias() = sigma_eff;
219
220 vars.temperature = N.dot(T);
222 {
223 auto const porosity =
225 .template value<double>(vars, x_position, t, dt);
226 auto const S_fr =
229 .template value<double>(vars, x_position, t, dt);
230 ip_data.phi_fr = S_fr * porosity;
231 }
232 }
233}

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::FrozenLiquid, MathLib::KelvinVector::Invariants< KelvinVectorSize >::identity2, NumLib::interpolateCoordinates(), localDOF(), MaterialPropertyLib::porosity, setInitialConditionsConcrete(), and MaterialPropertyLib::VariableArray::temperature.

Referenced by setInitialConditionsConcrete().

◆ setIPDataInitialConditions()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
std::size_t ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::setIPDataInitialConditions ( std::string_view const name,
double const * values,
int const integration_order )
overridevirtual

Returns number of read integration points.

Implements ProcessLib::ThermoHydroMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 107 of file ThermoHydroMechanicsFEM-impl.h.

110{
111 if (integration_order !=
112 static_cast<int>(_integration_method.getIntegrationOrder()))
113 {
114 OGS_FATAL(
115 "Setting integration point initial conditions; The integration "
116 "order of the local assembler for element {:d} is different from "
117 "the integration order in the initial condition.",
118 _element.getID());
119 }
120
121 if (name == "sigma")
122 {
123 if (_process_data.initial_stress.value)
124 {
125 OGS_FATAL(
126 "Setting initial conditions for stress from integration "
127 "point data and from a parameter '{:s}' is not possible "
128 "simultaneously.",
129 _process_data.initial_stress.value->name);
130 }
131
134 }
135 if (name == "epsilon_m")
136 {
139 }
140 if (name == "epsilon")
141 {
144 }
145 if (name.starts_with("material_state_variable_"))
146 {
147 name.remove_prefix(24);
148
149 // Using first ip data for solid material. TODO (naumov) move solid
150 // material into element, store only material state in IPs.
151 auto const& internal_variables =
152 _ip_data[0].solid_material.getInternalVariables();
153 if (auto const iv = std::find_if(
155 [&name](auto const& iv) { return iv.name == name; });
157 {
158 DBUG("Setting material state variable '{:s}'", name);
161 iv->reference);
162 }
163
164 int const element_id = _element.getID();
165 DBUG(
166 "The solid material of element {:d} (material ID {:d}) does not "
167 "have an internal state variable called {:s}.",
168 element_id, (*_process_data.material_ids)[element_id], name);
169 }
170
171 return 0;
172}
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:22
std::size_t setIntegrationPointDataMaterialStateVariables(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType member, std::function< std::span< double >(MaterialStateVariables &)> get_values_span)
std::size_t setIntegrationPointKelvinVectorData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)

References _element, _integration_method, _ip_data, _process_data, DBUG(), ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >::eps, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >::eps_m, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >::material_state_variables, OGS_FATAL, ProcessLib::setIntegrationPointDataMaterialStateVariables(), ProcessLib::setIntegrationPointKelvinVectorData(), and ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >::sigma_eff.

◆ setSigma()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
std::size_t ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::setSigma ( double const * values)
inlineprivate

◆ updateConstitutiveRelations()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
ConstitutiveRelationsValues< DisplacementDim > ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::updateConstitutiveRelations ( Eigen::Ref< Eigen::VectorXd const > const local_x,
Eigen::Ref< Eigen::VectorXd const > const local_x_prev,
ParameterLib::SpatialPosition const & x_position,
double const t,
double const dt,
IpData & ip_data,
IntegrationPointDataForOutput< DisplacementDim > & ip_data_output ) const
private

Definition at line 238 of file ThermoHydroMechanicsFEM-impl.h.

245{
246 assert(local_x.size() ==
248
249 auto const [T, p, u] = localDOF(local_x);
250 auto const [T_prev, p_prev, u_prev] = localDOF(local_x_prev);
251
252 auto const& solid_material =
254 _process_data.solid_materials, _process_data.material_ids,
255 _element.getID());
256
257 auto const& medium = _process_data.media_map.getMedium(_element.getID());
258 auto const& liquid_phase =
260 auto const& solid_phase =
262 auto* const frozen_liquid_phase =
265 : nullptr;
267
268 auto const& N_u = ip_data.N_u;
269 auto const& dNdx_u = ip_data.dNdx_u;
270
271 auto const& N = ip_data.N;
272 auto const& dNdx = ip_data.dNdx;
273
274 auto const T_int_pt = N.dot(T);
275 auto const T_prev_int_pt = N.dot(T_prev);
276 double const dT_int_pt = T_int_pt - T_prev_int_pt;
277
278 auto const x_coord =
279 x_position.getCoordinates().value()[0]; // r for axisymmetry
280 auto const B =
285
287
288 auto& eps = ip_data.eps;
289 eps.noalias() = B * u;
291 B * u_prev;
292
293 vars.temperature = T_int_pt;
294 double const p_int_pt = N.dot(p);
295 vars.liquid_phase_pressure = p_int_pt;
296 double const p_prev_int_pt = N.dot(p_prev);
297 double const dp_int_pt = p_int_pt - p_prev_int_pt;
298 vars.liquid_saturation = 1.0;
299
300 auto const solid_density =
302 .template value<double>(vars, x_position, t, dt);
303
304 auto const drho_SR_dT =
306 .template dValue<double>(vars,
308 x_position, t, dt);
309
310 auto const porosity =
312 .template value<double>(vars, x_position, t, dt);
313 vars.porosity = porosity;
314 ip_data.porosity = porosity;
315
316 crv.alpha_biot =
318 .template value<double>(vars, x_position, t, dt);
319 auto const& alpha = crv.alpha_biot;
320
321 auto const C_el = ip_data.computeElasticTangentStiffness(
322 t, x_position, dt, static_cast<double>(T_int_pt));
324 1 / solid_material.getBulkModulus(t, x_position, &C_el);
325
326 crv.beta_SR = (1 - alpha) * solid_skeleton_compressibility;
327
328 // Set mechanical variables for the intrinsic permeability model
329 // For stress dependent permeability.
330 {
331 auto const& identity2 = Invariants::identity2;
332 auto const sigma_total =
333 (ip_data.sigma_eff - alpha * p_int_pt * identity2).eval();
334 vars.total_stress.emplace<SymmetricTensor>(
336 }
337 // For strain dependent permeability
338 vars.volumetric_strain = Invariants::trace(ip_data.eps);
339 vars.equivalent_plastic_strain =
340 ip_data.material_state_variables->getEquivalentPlasticStrain();
341
342 auto const intrinsic_permeability =
345 .value(vars, x_position, t, dt));
346
347 auto const fluid_density =
349 .template value<double>(vars, x_position, t, dt);
350 ip_data_output.fluid_density = fluid_density;
351 vars.density = fluid_density;
352
353 auto const drho_dp =
355 .template dValue<double>(
357 x_position, t, dt);
358 crv.drho_LR_dp = drho_dp;
359
360 crv.fluid_compressibility = 1 / fluid_density * drho_dp;
361
362 crv.drho_LR_dT =
364 .template dValue<double>(vars,
366 x_position, t, dt);
367
371
372 // Use the viscosity model to compute the viscosity
373 ip_data_output.viscosity =
375 .template value<double>(vars, x_position, t, dt);
376 crv.K_over_mu = intrinsic_permeability / ip_data_output.viscosity;
377
378 crv.k_rel = 1.;
379 crv.dk_rel_dT = 0;
381 {
382 auto const S_fr =
383 (*medium)
385 .template value<double>(vars, x_position, t, dt);
386 ip_data.phi_fr = S_fr * porosity;
387
388 // The Sigmoid returns the frozen liquid (ice) saturation, i.e. the
389 // fraction of the pore space occupied by ice. dS_fr/dT is its
390 // temperature derivative; the absolute frozen volume fraction and its
391 // derivative are obtained by multiplying with the porosity.
392 auto const dS_fr_dT =
393 (*medium)
395 .template dValue<double>(
397 x_position, t, dt);
398
399 // Set frozen_liquid_saturation variable for the relative permeability
400 // calculation (fraction of pore space occupied by ice).
401 vars.frozen_liquid_saturation = S_fr;
402
403 crv.k_rel =
405 .property(
407 .template value<double>(vars, x_position, t, dt);
408
409 // dk_rel/dT = (dk_rel/dS_fr) * (dS_fr/dT)
410 auto const dk_rel_dS_fr =
412 .property(
414 .template dValue<double>(
415 vars,
417 x_position, t, dt);
418 crv.dk_rel_dT = dk_rel_dS_fr * dS_fr_dT;
419 }
420
421 auto const& b = _process_data.specific_body_force;
422
423 // Consider also anisotropic thermal expansion.
424 crv.solid_linear_thermal_expansion_coefficient =
427 .property(
429 .value(vars, x_position, t, dt));
430
433 crv.solid_linear_thermal_expansion_coefficient * dT_int_pt;
434
435 crv.K_pT_thermal_osmosis =
438 ip_data_output.viscosity);
439
441 -crv.k_rel * crv.K_over_mu * dNdx * p -
442 crv.K_pT_thermal_osmosis * dNdx * T +
443 (fluid_density * crv.k_rel) * crv.K_over_mu * b;
444 ip_data_output.velocity = velocity;
445 crv.dvelocity_dT =
446 -crv.dk_rel_dT * crv.K_over_mu * dNdx * p +
447 // TODO(naumov): - crv.K_pT_thermal_osmosis * dNdx * dT_dT +
448 (crv.dk_rel_dT * fluid_density + crv.k_rel * crv.drho_LR_dT) *
449 crv.K_over_mu * b;
450
451 //
452 // displacement equation, displacement part
453 //
454 auto& eps_m = ip_data.eps_m;
455 auto& eps_m_prev = ip_data.eps_m_prev;
456 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
457
458 crv.eps_v_dot = (vars.volumetric_strain - Invariants::trace(eps_prev)) / dt;
459
460 vars.mechanical_strain
462 eps_m);
463
464 crv.C = ip_data.updateConstitutiveRelation(vars, t, x_position, dt,
466
468
469 crv.beta =
471 (alpha - porosity) *
472 Invariants::trace(crv.solid_linear_thermal_expansion_coefficient);
473
474 //
475 // pressure equation, displacement part.
476 //
477 // Reusing Kup.transpose().
478
479 //
480 // temperature equation, temperature part.
481 //
482 crv.c_f =
485 .template value<double>(vars, x_position, t, dt);
486 crv.effective_thermal_conductivity =
488 medium
489 ->property(
491 .value(vars, x_position, t, dt));
492
494 medium
497 x_position, t, dt));
498
499 // Thermal conductivity is moved outside and zero matrix is passed instead
500 // due to multiplication with fluid's density times specific heat capacity.
501 crv.effective_thermal_conductivity.noalias() +=
502 fluid_density * crv.c_f *
504 _process_data.stabilizer, _element.getID(),
506 velocity, 0. /* phi */, 0. /* dispersivity_transversal */,
507 0. /*dispersivity_longitudinal*/);
508
509 double const c_s =
512 .template value<double>(vars, x_position, t, dt);
513
514 // Sensible part (no latent term here). Also modified by freezing terms.
515 crv.sensible_volumetric_heat_capacity =
516 porosity * fluid_density * crv.c_f +
517 (1.0 - porosity) * solid_density * c_s;
518 double dC_eff_dT = porosity * crv.drho_LR_dT * crv.c_f +
519 (1.0 - porosity) * drho_SR_dT * c_s;
520
522 {
524 double const phi_fr = ip_data.phi_fr;
525
526 auto const frozen_liquid_value =
528 {
529 return (*frozen_liquid_phase)[p].template value<double>(
530 vars, x_position, t, dt);
531 };
532
533 double const c_fr = frozen_liquid_value(
535
536 double const l_fr = frozen_liquid_value(
538
539 auto const dS_fr_dT =
540 (*medium)
542 .template dValue<double>(
544 x_position, t, dt);
545 double const dphi_fr_dT = dS_fr_dT * porosity;
546
547 auto const d2S_fr_dT2 =
548 (*medium)
550 .template d2Value<double>(
553 t, dt);
554 double const d2phi_fr_dT2 = d2S_fr_dT2 * porosity;
555
556 double const phi_fr_prev = [&]()
557 {
559 vars_prev.temperature = T_prev_int_pt;
560 auto const S_fr_prev =
563 .template value<double>(vars_prev, x_position, t, dt);
564 return S_fr_prev * porosity;
565 }();
566 ip_data.phi_fr_prev = phi_fr_prev;
567
568 double const rho_fr =
570 ip_data_output.rho_fr = rho_fr;
571
572 crv.rho += ip_data.phi_fr * rho_fr - ip_data.phi_fr * fluid_density;
573 crv.mass_exchange =
575 double const dmass_exchange_dT =
577 dphi_fr_dT * porosity * rho_fr * crv.drho_LR_dT /
579
580 // alpha_T^I
585 ->property(
587 .value(vars, x_position, t, dt));
588
592
593 crv.beta_T_SI =
594 porosity *
596 (alpha - porosity) *
598 crv.solid_linear_thermal_expansion_coefficient);
599
600 // alpha_{phi_I} -- linear expansion coeff. due to water-to-ice
601 // transition (phase change), and related phase_change_strain term
608 .value(vars, x_position, t, dt));
609
613
614 // eps0 ia a 'history variable' -- a solid matrix strain accrued
615 // prior to the onset of ice forming
616 auto& eps0 = ip_data.eps0;
617 auto const& eps0_prev = ip_data.eps0_prev;
618
619 // definition of eps_m_ice
620 auto& eps_m_ice = ip_data.eps_m_ice;
621 auto const& eps_m_ice_prev = ip_data.eps_m_ice_prev;
622
623 eps_m_ice.noalias() = eps_m_ice_prev + eps - eps_prev -
626
627 vars_ice.mechanical_strain
629 eps_m_ice);
630 auto const C_IR = ip_data.updateConstitutiveRelationIce(
631 *_process_data.ice_constitutive_relation, vars_ice, t, x_position,
633
634 auto const C_el_ice = ip_data.computeElasticTangentStiffnessIce(
635 *_process_data.ice_constitutive_relation, t, x_position, dt,
636 static_cast<double>(T_int_pt));
637 crv.beta_IR =
638 1. / _process_data.ice_constitutive_relation->getBulkModulus(
640
641 // Latent contribution L (effective = sensible - L).
642 crv.latent_volumetric_heat_capacity = l_fr * rho_fr * dphi_fr_dT;
643
644 // Freezing modifies only the sensible part; latent kept separate.
645 crv.sensible_volumetric_heat_capacity +=
646 -phi_fr * fluid_density * crv.c_f + phi_fr * rho_fr * c_fr;
647
648 crv.J_uu_fr = phi_fr * C_IR;
649
650 auto const& sigma_eff_ice = ip_data.sigma_eff_ice;
651 crv.r_u_fr = phi_fr * sigma_eff_ice;
652
654
655 // part of dMTT_dT derivative for freezing
657 phi_fr * crv.drho_LR_dT * crv.c_f +
659 double const storage_p_fr_coeff =
660 (porosity * crv.beta_IR + (alpha - porosity) * crv.beta_SR) *
662 (porosity * crv.fluid_compressibility +
663 (alpha - porosity) * crv.beta_SR);
664 crv.storage_p_fr = phi_fr / porosity * storage_p_fr_coeff;
665
666 double const dstorage_p_fr_coeff_dT =
667 (porosity * crv.beta_IR + (alpha - porosity) * crv.beta_SR) *
668 crv.drho_LR_dT / (fluid_density * fluid_density);
669
670 crv.J_pT_fr = (dphi_fr_dT * storage_p_fr_coeff +
673
674 crv.storage_T_fr =
675 phi_fr / porosity *
676 (crv.beta_T_SI * rho_fr / fluid_density - crv.beta) -
677 crv.mass_exchange;
678 double const dstorage_T_fr_dT =
680 (crv.beta_T_SI * rho_fr / fluid_density - crv.beta) +
681 phi_fr / porosity * crv.beta_T_SI * rho_fr * crv.drho_LR_dT /
684 crv.J_pT_fr += dstorage_T_fr_dT * dT_int_pt / dt;
685 }
686 crv.J_TT = dC_eff_dT * dT_int_pt / dt;
687 return crv;
688}
constexpr Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
MathLib::KelvinVector::KelvinVectorType< GlobalDim > formKelvinVector(MaterialPropertyLib::PropertyDataType const &values)
A function to form a Kelvin vector from strain or stress alike property like thermal expansivity for ...
double getLiquidThermalExpansivity(Phase const &phase, VariableArray const &vars, const double density, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
Eigen::MatrixXd computeHydrodynamicDispersion(NumericalStabilization const &stabilizer, std::size_t const element_id, Eigen::MatrixXd const &pore_diffusion_coefficient, Eigen::VectorXd const &velocity, double const porosity, double const solute_dispersivity_transverse, double const solute_dispersivity_longitudinal)
Eigen::Matrix< double, GlobalDim, GlobalDim > getThermoOsmoticCoefficient(MaterialPropertyLib::Medium const &medium, MaterialPropertyLib::VariableArray const &variable_array, ParameterLib::SpatialPosition const &pos, double const t, double const dt, Eigen::Matrix< double, GlobalDim, GlobalDim > const &intrinsic_permeability, double const liquid_dynamic_viscosity)
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.

References _element, _is_axially_symmetric, _process_data, MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::biot_coefficient, ProcessLib::LinearBMatrix::computeBMatrix(), ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::computeElasticTangentStiffness(), ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::computeElasticTangentStiffnessIce(), NumLib::computeHydrodynamicDispersion(), MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, displacement_size, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::dNdx, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::dNdx_u, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::eps, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::eps0, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::eps0_prev, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::eps_m, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::eps_m_ice, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::eps_m_ice_prev, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::eps_m_prev, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, ProcessLib::ThermoHydroMechanics::IntegrationPointDataForOutput< DisplacementDim >::fluid_density, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::formKelvinVector(), MaterialPropertyLib::frozen_liquid_saturation, MaterialPropertyLib::VariableArray::frozen_liquid_saturation, MaterialPropertyLib::FrozenLiquid, ParameterLib::SpatialPosition::getCoordinates(), MaterialPropertyLib::getLiquidThermalExpansivity(), ProcessLib::getThermoOsmoticCoefficient(), MathLib::KelvinVector::Invariants< KelvinVectorSize >::identity2, MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, localDOF(), ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::material_state_variables, MaterialPropertyLib::VariableArray::mechanical_strain, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::N, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::N_u, MaterialPropertyLib::permeability, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::phi_fr, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::phi_fr_prev, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::porosity, pressure_size, MaterialPropertyLib::relative_permeability, ProcessLib::ThermoHydroMechanics::IntegrationPointDataForOutput< DisplacementDim >::rho_fr, MaterialLib::Solids::selectSolidConstitutiveRelation(), ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::sigma_eff, ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::sigma_eff_ice, MaterialPropertyLib::Solid, MaterialPropertyLib::specific_heat_capacity, MaterialPropertyLib::specific_latent_heat, MaterialPropertyLib::temperature, MaterialPropertyLib::VariableArray::temperature, temperature_size, MaterialPropertyLib::thermal_conductivity, MaterialPropertyLib::thermal_expansivity, MaterialPropertyLib::VariableArray::total_stress, MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace(), ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::updateConstitutiveRelation(), ProcessLib::ThermoHydroMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::updateConstitutiveRelationIce(), ProcessLib::ThermoHydroMechanics::IntegrationPointDataForOutput< DisplacementDim >::velocity, MaterialPropertyLib::viscosity, ProcessLib::ThermoHydroMechanics::IntegrationPointDataForOutput< DisplacementDim >::viscosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

Referenced by assembleWithJacobian(), and postTimestepConcrete().

Member Data Documentation

◆ _element

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
MeshLib::Element const& ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_element
private

◆ _integration_method

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
NumLib::GenericIntegrationMethod const& ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_integration_method
private

◆ _ip_data

◆ _ip_data_output

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
std::vector<IntegrationPointDataForOutput<DisplacementDim>, Eigen::aligned_allocator< IntegrationPointDataForOutput<DisplacementDim> > > ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data_output
private

◆ _is_axially_symmetric

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
bool const ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_is_axially_symmetric
private

◆ _process_data

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
ThermoHydroMechanicsProcessData<DisplacementDim>& ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_process_data
private

◆ _secondary_data

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType> ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_secondary_data
private

◆ displacement_index

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
const int ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::displacement_index = ShapeFunctionPressure::NPOINTS * 2
staticprivate

Definition at line 454 of file ThermoHydroMechanicsFEM.h.

Referenced by assembleWithJacobian().

◆ displacement_size

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
const int ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::displacement_size
staticprivate
Initial value:
=
ShapeFunctionDisplacement::NPOINTS * DisplacementDim

Definition at line 455 of file ThermoHydroMechanicsFEM.h.

Referenced by assembleWithJacobian(), getNumberOfVectorElementsForDeformation(), and updateConstitutiveRelations().

◆ KelvinVectorSize

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
int const ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::KelvinVectorSize
static
Initial value:

Definition at line 62 of file ThermoHydroMechanicsFEM.h.

◆ N_u_op

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
auto& ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::N_u_op
staticconstexpr
Initial value:
DisplacementDim,
constexpr Eigen::CwiseNullaryOp< EigenBlockMatrixViewFunctor< D, M >, typename EigenBlockMatrixViewFunctor< D, M >::Matrix > eigenBlockMatrixView(const Eigen::MatrixBase< M > &matrix)
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType

Definition at line 68 of file ThermoHydroMechanicsFEM.h.

Referenced by assembleWithJacobian().

◆ pressure_index

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
const int ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::pressure_index = ShapeFunctionPressure::NPOINTS
staticprivate

◆ pressure_size

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
const int ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::pressure_size = ShapeFunctionPressure::NPOINTS
staticprivate

Definition at line 453 of file ThermoHydroMechanicsFEM.h.

Referenced by assembleWithJacobian(), and updateConstitutiveRelations().

◆ temperature_index

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
const int ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::temperature_index = 0
staticprivate

◆ temperature_size

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
const int ProcessLib::ThermoHydroMechanics::ThermoHydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::temperature_size = ShapeFunctionPressure::NPOINTS
staticprivate

Definition at line 451 of file ThermoHydroMechanicsFEM.h.

Referenced by assembleWithJacobian(), and updateConstitutiveRelations().


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