Loading web-font TeX/Main/Regular
OGS
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > Class Template Reference

Detailed Description

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
class ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >

Definition at line 51 of file RichardsMechanicsFEM.h.

#include <RichardsMechanicsFEM.h>

Inheritance diagram for ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:
[legend]
Collaboration diagram for ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:
[legend]

Public Types

using ShapeMatricesTypeDisplacement
 
using ShapeMatricesTypePressure
 
using GlobalDimMatrixType
 
using BMatricesType
 
using KelvinVectorType = typename BMatricesType::KelvinVectorType
 
using IpData
 
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>
 
using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>
 

Public Member Functions

 RichardsMechanicsLocalAssembler (RichardsMechanicsLocalAssembler const &)=delete
 
 RichardsMechanicsLocalAssembler (RichardsMechanicsLocalAssembler &&)=delete
 
 RichardsMechanicsLocalAssembler (MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, RichardsMechanicsProcessData< DisplacementDim > &process_data)
 
void setInitialConditionsConcrete (Eigen::VectorXd const local_x, double const t, int const process_id) override
 
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_rhs_data) 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 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) override
 
void initializeConcrete () 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.
 
- Public Member Functions inherited from ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >
 LocalAssemblerInterface (MeshLib::Element const &e, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, RichardsMechanicsProcessData< DisplacementDim > &process_data)
 
std::size_t setIPDataInitialConditions (std::string_view name, double const *values, int const integration_order)
 
std::vector< double > getMaterialStateVariableInternalState (std::function< std::span< double >(typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables &)> const &get_values_span, int const &n_components) const
 
unsigned getNumberOfIntegrationPoints () const
 
int getMaterialID () const
 
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt (unsigned integration_point) const
 
void postTimestepConcrete (Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override final
 
- 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 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 Member Functions

void assembleWithJacobianForDeformationEquations (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
 
void assembleWithJacobianForPressureEquations (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
 

Static Private Member Functions

static void assembleWithJacobianEvalConstitutiveSetting (double const t, double const dt, ParameterLib::SpatialPosition const &x_position, IpData &ip_data, MPL::VariableArray &variables, MPL::VariableArray &variables_prev, MPL::Medium const *const medium, TemperatureData const T_data, CapillaryPressureData< DisplacementDim > const &p_cap_data, ConstitutiveData< DisplacementDim > &CD, StatefulData< DisplacementDim > &SD, StatefulDataPrev< DisplacementDim > const &SD_prev, std::optional< MicroPorosityParameters > const &micro_porosity_parameters, MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material, ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > &material_state_data)
 

Private Attributes

std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
 
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
 

Static Private Attributes

static const int pressure_index = 0
 
static const int pressure_size = ShapeFunctionPressure::NPOINTS
 
static const int displacement_index = ShapeFunctionPressure::NPOINTS
 
static const int displacement_size
 

Additional Inherited Members

- Static Public Member Functions inherited from ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >
static auto getReflectionDataForOutput ()
 
- Protected Attributes inherited from ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >
RichardsMechanicsProcessData< DisplacementDim > & process_data_
 
NumLib::GenericIntegrationMethod const & integration_method_
 
MeshLib::Element const & element_
 
bool const is_axially_symmetric_
 
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material_
 
std::vector< StatefulData< DisplacementDim > > current_states_
 
std::vector< StatefulDataPrev< DisplacementDim > > prev_states_
 
std::vector< ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > > material_states_
 
std::vector< OutputData< DisplacementDim > > output_data_
 

Member Typedef Documentation

◆ BMatricesType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::BMatricesType
Initial value:
BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 63 of file RichardsMechanicsFEM.h.

◆ GlobalDimMatrixType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::GlobalDimMatrixType
Initial value:

Definition at line 60 of file RichardsMechanicsFEM.h.

◆ Invariants

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

Definition at line 74 of file RichardsMechanicsFEM.h.

◆ IpData

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

Definition at line 67 of file RichardsMechanicsFEM.h.

◆ KelvinVectorType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::KelvinVectorType = typename BMatricesType::KelvinVectorType

Definition at line 65 of file RichardsMechanicsFEM.h.

◆ ShapeMatricesTypeDisplacement

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypeDisplacement

◆ ShapeMatricesTypePressure

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypePressure

◆ SymmetricTensor

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

Definition at line 76 of file RichardsMechanicsFEM.h.

Constructor & Destructor Documentation

◆ RichardsMechanicsLocalAssembler() [1/3]

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

◆ RichardsMechanicsLocalAssembler() [2/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::RichardsMechanicsLocalAssembler ( RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > && )
delete

◆ RichardsMechanicsLocalAssembler() [3/3]

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

Definition at line 122 of file RichardsMechanicsFEM-impl.h.

130 e, integration_method, is_axially_symmetric, process_data}
131{
132 unsigned const n_integration_points =
134
135 _ip_data.resize(n_integration_points);
136 _secondary_data.N_u.resize(n_integration_points);
137
138 auto const shape_matrices_u =
139 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
141 DisplacementDim>(e, is_axially_symmetric,
142 this->integration_method_);
143
144 auto const shape_matrices_p =
145 NumLib::initShapeMatrices<ShapeFunctionPressure,
146 ShapeMatricesTypePressure, DisplacementDim>(
147 e, is_axially_symmetric, this->integration_method_);
148
149 auto const& medium =
150 this->process_data_.media_map.getMedium(this->element_.getID());
151
153 x_position.setElementID(this->element_.getID());
154 for (unsigned ip = 0; ip < n_integration_points; ip++)
155 {
156 x_position.setIntegrationPoint(ip);
157 auto& ip_data = _ip_data[ip];
158 auto const& sm_u = shape_matrices_u[ip];
159 _ip_data[ip].integration_weight =
161 sm_u.integralMeasure * sm_u.detJ;
162
163 ip_data.N_u = sm_u.N;
164 ip_data.dNdx_u = sm_u.dNdx;
165
166 ip_data.N_p = shape_matrices_p[ip].N;
167 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
168
169 // Initial porosity. Could be read from integration point data or mesh.
170 auto& porosity =
171 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
172 this->current_states_[ip])
173 .phi;
174 porosity = medium->property(MPL::porosity)
175 .template initialValue<double>(
176 x_position,
177 std::numeric_limits<
178 double>::quiet_NaN() /* t independent */);
179
180 auto& transport_porosity =
181 std::get<
183 this->current_states_[ip])
184 .phi;
186 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
187 {
189 medium->property(MPL::transport_porosity)
190 .template initialValue<double>(
191 x_position,
192 std::numeric_limits<
193 double>::quiet_NaN() /* t independent */);
194 }
195
196 _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
197 }
198}
double getWeight() const
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::vector< StatefulData< DisplacementDim > > current_states_
LocalAssemblerInterface(MeshLib::Element const &e, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, RichardsMechanicsProcessData< DisplacementDim > &process_data)
RichardsMechanicsProcessData< DisplacementDim > & process_data_
NumLib::GenericIntegrationMethod const & integration_method_
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u

Member Function Documentation

◆ assemble()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::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_rhs_data )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 342 of file RichardsMechanicsFEM-impl.h.

348{
349 assert(local_x.size() == pressure_size + displacement_size);
350
351 auto p_L =
352 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
353 pressure_size> const>(local_x.data() + pressure_index,
355
356 auto u =
357 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
358 displacement_size> const>(local_x.data() + displacement_index,
360
361 auto p_L_prev =
362 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
363 pressure_size> const>(local_x_prev.data() + pressure_index,
365
366 auto u_prev =
367 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
368 displacement_size> const>(local_x_prev.data() + displacement_index,
370
372 typename ShapeMatricesTypeDisplacement::template MatrixType<
375 local_K_data, displacement_size + pressure_size,
377
379 typename ShapeMatricesTypeDisplacement::template MatrixType<
382 local_M_data, displacement_size + pressure_size,
384
386 typename ShapeMatricesTypeDisplacement::template VectorType<
388 local_rhs_data, displacement_size + pressure_size);
389
390 auto const& identity2 = MathLib::KelvinVector::Invariants<
392 DisplacementDim)>::identity2;
393
394 auto const& medium =
395 this->process_data_.media_map.getMedium(this->element_.getID());
396 auto const& liquid_phase = medium->phase("AqueousLiquid");
397 auto const& solid_phase = medium->phase("Solid");
398 MPL::VariableArray variables;
399 MPL::VariableArray variables_prev;
400
402 x_position.setElementID(this->element_.getID());
403
404 unsigned const n_integration_points =
406 for (unsigned ip = 0; ip < n_integration_points; ip++)
407 {
408 x_position.setIntegrationPoint(ip);
409 auto const& w = _ip_data[ip].integration_weight;
410
411 auto const& N_u = _ip_data[ip].N_u;
412 auto const& dNdx_u = _ip_data[ip].dNdx_u;
413
414 auto const& N_p = _ip_data[ip].N_p;
415 auto const& dNdx_p = _ip_data[ip].dNdx_p;
416
417 auto const x_coord =
418 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
420 this->element_, N_u);
421 auto const B =
422 LinearBMatrix::computeBMatrix<DisplacementDim,
423 ShapeFunctionDisplacement::NPOINTS,
425 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
426
427 auto& eps =
428 std::get<StrainData<DisplacementDim>>(this->current_states_[ip]);
429 eps.eps.noalias() = B * u;
430
431 auto& S_L =
432 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
433 this->current_states_[ip])
434 .S_L;
435 auto const S_L_prev =
436 std::get<
437 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
438 this->prev_states_[ip])
439 ->S_L;
440
441 double p_cap_ip;
442 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
443
444 double p_cap_prev_ip;
445 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
446
447 variables.capillary_pressure = p_cap_ip;
448 variables.liquid_phase_pressure = -p_cap_ip;
449 // setting pG to 1 atm
450 // TODO : rewrite equations s.t. p_L = pG-p_cap
451 variables.gas_phase_pressure = 1.0e5;
452
453 auto const temperature =
454 medium->property(MPL::PropertyType::reference_temperature)
455 .template value<double>(variables, x_position, t, dt);
456 variables.temperature = temperature;
457
458 auto const alpha =
459 medium->property(MPL::PropertyType::biot_coefficient)
460 .template value<double>(variables, x_position, t, dt);
461 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
462 t, x_position, dt, temperature, this->solid_material_,
463 *this->material_states_[ip].material_state_variables);
464
465 auto const beta_SR = (1 - alpha) / this->solid_material_.getBulkModulus(
466 t, x_position, &C_el);
467 variables.grain_compressibility = beta_SR;
468
469 auto const rho_LR =
470 liquid_phase.property(MPL::PropertyType::density)
471 .template value<double>(variables, x_position, t, dt);
472 variables.density = rho_LR;
473
474 auto const& b = this->process_data_.specific_body_force;
475
476 S_L = medium->property(MPL::PropertyType::saturation)
477 .template value<double>(variables, x_position, t, dt);
478 variables.liquid_saturation = S_L;
479 variables_prev.liquid_saturation = S_L_prev;
480
481 // tangent derivative for Jacobian
482 double const dS_L_dp_cap =
483 medium->property(MPL::PropertyType::saturation)
484 .template dValue<double>(variables,
485 MPL::Variable::capillary_pressure,
486 x_position, t, dt);
487 // secant derivative from time discretization for storage
488 // use tangent, if secant is not available
489 double const DeltaS_L_Deltap_cap =
490 (p_cap_ip == p_cap_prev_ip)
491 ? dS_L_dp_cap
492 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
493
494 auto const chi = [medium, x_position, t, dt](double const S_L)
495 {
497 vs.liquid_saturation = S_L;
498 return medium->property(MPL::PropertyType::bishops_effective_stress)
499 .template value<double>(vs, x_position, t, dt);
500 };
501 double const chi_S_L = chi(S_L);
502 double const chi_S_L_prev = chi(S_L_prev);
503
504 double const p_FR = -chi_S_L * p_cap_ip;
505 variables.effective_pore_pressure = p_FR;
506 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
507
508 // Set volumetric strain rate for the general case without swelling.
509 variables.volumetric_strain = Invariants::trace(eps.eps);
510 variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
511
512 auto& phi = std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
513 this->current_states_[ip])
514 .phi;
515 { // Porosity update
516 auto const phi_prev = std::get<PrevState<
518 this->prev_states_[ip])
519 ->phi;
520 variables_prev.porosity = phi_prev;
521 phi = medium->property(MPL::PropertyType::porosity)
522 .template value<double>(variables, variables_prev,
523 x_position, t, dt);
524 variables.porosity = phi;
525 }
526
527 if (alpha < phi)
528 {
529 OGS_FATAL(
530 "RichardsMechanics: Biot-coefficient {} is smaller than "
531 "porosity {} in element/integration point {}/{}.",
532 alpha, phi, this->element_.getID(), ip);
533 }
534
535 // Swelling and possibly volumetric strain rate update.
536 {
537 auto& sigma_sw =
538 std::get<ProcessLib::ThermoRichardsMechanics::
539 ConstitutiveStress_StrainTemperature::
540 SwellingDataStateful<DisplacementDim>>(
541 this->current_states_[ip])
542 .sigma_sw;
543 auto const& sigma_sw_prev = std::get<PrevState<
544 ProcessLib::ThermoRichardsMechanics::
545 ConstitutiveStress_StrainTemperature::SwellingDataStateful<
546 DisplacementDim>>>(this->prev_states_[ip])
547 ->sigma_sw;
548
549 // If there is swelling, compute it. Update volumetric strain rate,
550 // s.t. it corresponds to the mechanical part only.
551 sigma_sw = sigma_sw_prev;
552 if (solid_phase.hasProperty(
553 MPL::PropertyType::swelling_stress_rate))
554 {
555 auto const sigma_sw_dot =
558 solid_phase[MPL::PropertyType::swelling_stress_rate]
559 .value(variables, variables_prev, x_position, t,
560 dt)));
561 sigma_sw += sigma_sw_dot * dt;
562
563 // !!! Misusing volumetric strain for mechanical volumetric
564 // strain just to update the transport porosity !!!
565 variables.volumetric_strain +=
566 identity2.transpose() * C_el.inverse() * sigma_sw;
567 variables_prev.volumetric_strain +=
568 identity2.transpose() * C_el.inverse() * sigma_sw_prev;
569 }
570
571 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
572 {
573 auto& transport_porosity =
574 std::get<ProcessLib::ThermoRichardsMechanics::
575 TransportPorosityData>(
576 this->current_states_[ip])
577 .phi;
578 auto const transport_porosity_prev =
579 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
580 TransportPorosityData>>(
581 this->prev_states_[ip])
582 ->phi;
583 variables_prev.transport_porosity = transport_porosity_prev;
584
586 medium->property(MPL::PropertyType::transport_porosity)
587 .template value<double>(variables, variables_prev,
588 x_position, t, dt);
590 }
591 else
592 {
593 variables.transport_porosity = phi;
594 }
595 }
596
597 double const k_rel =
598 medium->property(MPL::PropertyType::relative_permeability)
599 .template value<double>(variables, x_position, t, dt);
600 auto const mu =
601 liquid_phase.property(MPL::PropertyType::viscosity)
602 .template value<double>(variables, x_position, t, dt);
603
604 auto const& sigma_sw =
605 std::get<ProcessLib::ThermoRichardsMechanics::
606 ConstitutiveStress_StrainTemperature::
607 SwellingDataStateful<DisplacementDim>>(
608 this->current_states_[ip])
609 .sigma_sw;
610 auto const& sigma_eff =
611 std::get<ProcessLib::ThermoRichardsMechanics::
612 ConstitutiveStress_StrainTemperature::
613 EffectiveStressData<DisplacementDim>>(
614 this->current_states_[ip])
615 .sigma_eff;
616
617 // Set mechanical variables for the intrinsic permeability model
618 // For stress dependent permeability.
619 {
620 auto const sigma_total =
621 (sigma_eff - alpha * p_FR * identity2).eval();
622
623 // For stress dependent permeability.
624 variables.total_stress.emplace<SymmetricTensor>(
626 sigma_total));
627 }
628
629 variables.equivalent_plastic_strain =
630 this->material_states_[ip]
631 .material_state_variables->getEquivalentPlasticStrain();
632
633 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
634 medium->property(MPL::PropertyType::permeability)
635 .value(variables, x_position, t, dt));
636
637 GlobalDimMatrixType const rho_K_over_mu =
638 K_intrinsic * rho_LR * k_rel / mu;
639
640 //
641 // displacement equation, displacement part
642 //
643 {
644 auto& eps_m =
645 std::get<ProcessLib::ThermoRichardsMechanics::
646 ConstitutiveStress_StrainTemperature::
647 MechanicalStrainData<DisplacementDim>>(
648 this->current_states_[ip])
649 .eps_m;
650 eps_m.noalias() =
651 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
652 ? eps.eps + C_el.inverse() * sigma_sw
653 : eps.eps;
654 variables.mechanical_strain.emplace<
656 eps_m);
657 }
658
659 {
660 auto& SD = this->current_states_[ip];
661 auto const& SD_prev = this->prev_states_[ip];
662 auto& sigma_eff =
663 std::get<ProcessLib::ThermoRichardsMechanics::
664 ConstitutiveStress_StrainTemperature::
665 EffectiveStressData<DisplacementDim>>(SD);
666 auto const& sigma_eff_prev = std::get<
667 PrevState<ProcessLib::ThermoRichardsMechanics::
668 ConstitutiveStress_StrainTemperature::
669 EffectiveStressData<DisplacementDim>>>(
670 SD_prev);
671 auto const& eps_m =
672 std::get<ProcessLib::ThermoRichardsMechanics::
673 ConstitutiveStress_StrainTemperature::
674 MechanicalStrainData<DisplacementDim>>(SD);
675 auto& eps_m_prev = std::get<
676 PrevState<ProcessLib::ThermoRichardsMechanics::
677 ConstitutiveStress_StrainTemperature::
678 MechanicalStrainData<DisplacementDim>>>(
679 SD_prev);
680
681 _ip_data[ip].updateConstitutiveRelation(
682 variables, t, x_position, dt, temperature, sigma_eff,
683 sigma_eff_prev, eps_m, eps_m_prev, this->solid_material_,
684 this->material_states_[ip].material_state_variables);
685 }
686
687 // p_SR
688 variables.solid_grain_pressure =
689 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
690 auto const rho_SR =
691 solid_phase.property(MPL::PropertyType::density)
692 .template value<double>(variables, x_position, t, dt);
693
694 //
695 // displacement equation, displacement part
696 //
697 double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
698 rhs.template segment<displacement_size>(displacement_index).noalias() -=
699 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
700
701 //
702 // pressure equation, pressure part.
703 //
704 auto const beta_LR =
705 1 / rho_LR *
706 liquid_phase.property(MPL::PropertyType::density)
707 .template dValue<double>(variables,
708 MPL::Variable::liquid_phase_pressure,
709 x_position, t, dt);
710
711 double const a0 = S_L * (alpha - phi) * beta_SR;
712 // Volumetric average specific storage of the solid and fluid phases.
713 double const specific_storage =
714 DeltaS_L_Deltap_cap * (p_cap_ip * a0 - phi) +
715 S_L * (phi * beta_LR + a0);
716 M.template block<pressure_size, pressure_size>(pressure_index,
718 .noalias() += N_p.transpose() * rho_LR * specific_storage * N_p * w;
719
720 K.template block<pressure_size, pressure_size>(pressure_index,
722 .noalias() += dNdx_p.transpose() * rho_K_over_mu * dNdx_p * w;
723
724 rhs.template segment<pressure_size>(pressure_index).noalias() +=
725 dNdx_p.transpose() * rho_LR * rho_K_over_mu * b * w;
726
727 //
728 // displacement equation, pressure part
729 //
730 K.template block<displacement_size, pressure_size>(displacement_index,
732 .noalias() -= B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
733
734 //
735 // pressure equation, displacement part.
736 //
737 M.template block<pressure_size, displacement_size>(pressure_index,
739 .noalias() += N_p.transpose() * S_L * rho_LR * alpha *
740 identity2.transpose() * B * w;
741 }
742
743 if (this->process_data_.apply_mass_lumping)
744 {
745 auto Mpp = M.template block<pressure_size, pressure_size>(
747 Mpp = Mpp.colwise().sum().eval().asDiagonal();
748 }
749}
#define OGS_FATAL(...)
Definition Error.h:26
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
@ rho
density. For some models, rho substitutes p
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
KelvinVectorType< DisplacementDim > tensorToKelvin(Eigen::Matrix< double, 3, 3 > const &m)
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
auto eval(Function &f, Tuples &... ts) -> typename detail::GetFunctionReturnType< decltype(&Function::eval)>::type
Definition Apply.h:274
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.
std::vector< StatefulDataPrev< DisplacementDim > > prev_states_
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material_
std::vector< ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > > material_states_

References MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::formEigenTensor< 3 >(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::VariableArray::grain_compressibility, NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::VariableArray::mechanical_strain, OGS_FATAL, MaterialPropertyLib::VariableArray::porosity, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::temperature, MathLib::KelvinVector::tensorToKelvin(), MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::VariableArray::transport_porosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ assembleWithJacobian()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< 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 1071 of file RichardsMechanicsFEM-impl.h.

1077{
1078 assert(local_x.size() == pressure_size + displacement_size);
1079
1080 auto p_L =
1081 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
1082 pressure_size> const>(local_x.data() + pressure_index,
1084
1085 auto u =
1086 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
1087 displacement_size> const>(local_x.data() + displacement_index,
1089
1090 auto p_L_prev =
1091 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
1092 pressure_size> const>(local_x_prev.data() + pressure_index,
1094 auto u_prev =
1095 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
1096 displacement_size> const>(local_x_prev.data() + displacement_index,
1098
1099 auto local_Jac = MathLib::createZeroedMatrix<
1100 typename ShapeMatricesTypeDisplacement::template MatrixType<
1103 local_Jac_data, displacement_size + pressure_size,
1105
1106 auto local_rhs = MathLib::createZeroedVector<
1107 typename ShapeMatricesTypeDisplacement::template VectorType<
1109 local_rhs_data, displacement_size + pressure_size);
1110
1111 auto const& identity2 = MathLib::KelvinVector::Invariants<
1113 DisplacementDim)>::identity2;
1114
1116 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1118
1119 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_p =
1120 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1122
1123 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S_Jpp =
1124 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1126
1127 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S =
1128 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1130
1131 typename ShapeMatricesTypeDisplacement::template MatrixType<
1133 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
1136
1137 typename ShapeMatricesTypeDisplacement::template MatrixType<
1139 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
1142
1143 auto const& medium =
1144 this->process_data_.media_map.getMedium(this->element_.getID());
1145 auto const& liquid_phase = medium->phase("AqueousLiquid");
1146 auto const& solid_phase = medium->phase("Solid");
1147 MPL::VariableArray variables;
1148 MPL::VariableArray variables_prev;
1149
1151 x_position.setElementID(this->element_.getID());
1152
1153 unsigned const n_integration_points =
1155 for (unsigned ip = 0; ip < n_integration_points; ip++)
1156 {
1158 auto& SD = this->current_states_[ip];
1159 auto const& SD_prev = this->prev_states_[ip];
1160 [[maybe_unused]] auto models = createConstitutiveModels(
1161 this->process_data_, this->solid_material_);
1162
1163 x_position.setIntegrationPoint(ip);
1164 auto const& w = _ip_data[ip].integration_weight;
1165
1166 auto const& N_u = _ip_data[ip].N_u;
1167 auto const& dNdx_u = _ip_data[ip].dNdx_u;
1168
1169 auto const& N_p = _ip_data[ip].N_p;
1170 auto const& dNdx_p = _ip_data[ip].dNdx_p;
1171
1172 auto const x_coord =
1173 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1175 this->element_, N_u);
1176 auto const B =
1177 LinearBMatrix::computeBMatrix<DisplacementDim,
1178 ShapeFunctionDisplacement::NPOINTS,
1180 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
1181
1182 double p_cap_ip;
1183 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
1184
1185 double p_cap_prev_ip;
1186 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
1187
1188 variables.capillary_pressure = p_cap_ip;
1189 variables.liquid_phase_pressure = -p_cap_ip;
1190 // setting pG to 1 atm
1191 // TODO : rewrite equations s.t. p_L = pG-p_cap
1192 variables.gas_phase_pressure = 1.0e5;
1193
1194 auto const temperature =
1195 medium->property(MPL::PropertyType::reference_temperature)
1196 .template value<double>(variables, x_position, t, dt);
1197 variables.temperature = temperature;
1198
1199 std::get<StrainData<DisplacementDim>>(SD).eps.noalias() = B * u;
1200
1202 t, dt, x_position, _ip_data[ip], variables, variables_prev, medium,
1204 CapillaryPressureData<DisplacementDim>{
1205 p_cap_ip, p_cap_prev_ip,
1206 Eigen::Vector<double, DisplacementDim>::Zero()},
1207 CD, SD, SD_prev, this->process_data_.micro_porosity_parameters,
1208 this->solid_material_, this->material_states_[ip]);
1209
1210 {
1211 auto const& C = *std::get<StiffnessTensor<DisplacementDim>>(CD);
1212 local_Jac
1213 .template block<displacement_size, displacement_size>(
1215 .noalias() += B.transpose() * C * B * w;
1216 }
1217
1218 auto const& b = this->process_data_.specific_body_force;
1219
1220 {
1221 auto const& sigma_eff =
1222 std::get<ProcessLib::ThermoRichardsMechanics::
1223 ConstitutiveStress_StrainTemperature::
1224 EffectiveStressData<DisplacementDim>>(
1225 this->current_states_[ip])
1226 .sigma_eff;
1227 double const rho = *std::get<Density>(CD);
1228 local_rhs.template segment<displacement_size>(displacement_index)
1229 .noalias() -= (B.transpose() * sigma_eff -
1230 N_u_op(N_u).transpose() * rho * b) *
1231 w;
1232 }
1233
1234 //
1235 // displacement equation, pressure part
1236 //
1237
1238 double const alpha =
1239 *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD);
1240 double const dS_L_dp_cap =
1241 std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(
1242 CD)
1243 .dS_L_dp_cap;
1244
1245 {
1246 double const chi_S_L =
1247 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1248 .chi_S_L;
1249 Kup.noalias() +=
1250 B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
1251 double const dchi_dS_L =
1252 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1253 .dchi_dS_L;
1254
1255 local_Jac
1256 .template block<displacement_size, pressure_size>(
1258 .noalias() -= B.transpose() * alpha *
1259 (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
1260 identity2 * N_p * w;
1261 }
1262
1263 double const phi =
1264 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi;
1265 double const rho_LR = *std::get<LiquidDensity>(CD);
1266 local_Jac
1267 .template block<displacement_size, pressure_size>(
1269 .noalias() +=
1270 N_u_op(N_u).transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
1271
1272 // For the swelling stress with double structure model the corresponding
1273 // Jacobian u-p entry would be required, but it does not improve
1274 // convergence and sometimes worsens it:
1275 // if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1276 // {
1277 // -B.transpose() *
1278 // dsigma_sw_dS_L_m* dS_L_m_dp_cap_m*(p_L_m - p_L_m_prev) /
1279 // (p_cap_ip - p_cap_prev_ip) * N_p* w;
1280 // }
1281 if (!medium->hasProperty(MPL::PropertyType::saturation_micro) &&
1282 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
1283 {
1284 using DimMatrix = Eigen::Matrix<double, 3, 3>;
1285 auto const dsigma_sw_dS_L =
1287 solid_phase
1288 .property(MPL::PropertyType::swelling_stress_rate)
1289 .template dValue<DimMatrix>(
1290 variables, variables_prev,
1291 MPL::Variable::liquid_saturation, x_position, t,
1292 dt));
1293 local_Jac
1294 .template block<displacement_size, pressure_size>(
1296 .noalias() +=
1297 B.transpose() * dsigma_sw_dS_L * dS_L_dp_cap * N_p * w;
1298 }
1299 //
1300 // pressure equation, displacement part.
1301 //
1302 double const S_L =
1303 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
1304 this->current_states_[ip])
1305 .S_L;
1306 if (this->process_data_.explicit_hm_coupling_in_unsaturated_zone)
1307 {
1308 double const chi_S_L_prev = std::get<PrevState<
1310 ->chi_S_L;
1311 Kpu.noalias() += N_p.transpose() * chi_S_L_prev * rho_LR * alpha *
1312 identity2.transpose() * B * w;
1313 }
1314 else
1315 {
1316 Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
1317 identity2.transpose() * B * w;
1318 }
1319
1320 //
1321 // pressure equation, pressure part.
1322 //
1323
1324 double const k_rel =
1326 DisplacementDim>>(CD)
1327 .k_rel;
1328 auto const& K_intrinsic =
1330 DisplacementDim>>(CD)
1331 .Ki;
1332 double const mu =
1333 *std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(
1334 CD);
1335
1336 GlobalDimMatrixType const rho_Ki_over_mu = K_intrinsic * rho_LR / mu;
1337
1338 laplace_p.noalias() +=
1339 dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
1340
1341 auto const beta_LR =
1342 1 / rho_LR *
1343 liquid_phase.property(MPL::PropertyType::density)
1344 .template dValue<double>(variables,
1345 MPL::Variable::liquid_phase_pressure,
1346 x_position, t, dt);
1347
1348 double const beta_SR =
1349 std::get<
1351 CD)
1352 .beta_SR;
1353 double const a0 = (alpha - phi) * beta_SR;
1354 double const specific_storage_a_p = S_L * (phi * beta_LR + S_L * a0);
1355 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
1356
1357 double const dspecific_storage_a_p_dp_cap =
1358 dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0);
1359 double const dspecific_storage_a_S_dp_cap =
1360 -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
1361
1362 storage_p_a_p.noalias() +=
1363 N_p.transpose() * rho_LR * specific_storage_a_p * N_p * w;
1364
1365 double const DeltaS_L_Deltap_cap =
1366 std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap;
1367 storage_p_a_S.noalias() -= N_p.transpose() * rho_LR *
1368 specific_storage_a_S * DeltaS_L_Deltap_cap *
1369 N_p * w;
1370
1371 local_Jac
1372 .template block<pressure_size, pressure_size>(pressure_index,
1374 .noalias() += N_p.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
1375 rho_LR * dspecific_storage_a_p_dp_cap * N_p * w;
1376
1377 double const S_L_prev =
1378 std::get<
1379 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
1380 this->prev_states_[ip])
1381 ->S_L;
1382 storage_p_a_S_Jpp.noalias() -=
1383 N_p.transpose() * rho_LR *
1384 ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
1385 specific_storage_a_S * dS_L_dp_cap) /
1386 dt * N_p * w;
1387
1388 if (!this->process_data_.explicit_hm_coupling_in_unsaturated_zone)
1389 {
1390 local_Jac
1391 .template block<pressure_size, pressure_size>(pressure_index,
1393 .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap * alpha *
1394 identity2.transpose() * B * (u - u_prev) / dt *
1395 N_p * w;
1396 }
1397
1398 double const dk_rel_dS_l =
1399 medium->property(MPL::PropertyType::relative_permeability)
1400 .template dValue<double>(variables,
1401 MPL::Variable::liquid_saturation,
1402 x_position, t, dt);
1404 grad_p_cap = -dNdx_p * p_L;
1405 local_Jac
1406 .template block<pressure_size, pressure_size>(pressure_index,
1408 .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
1409 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1410
1411 local_Jac
1412 .template block<pressure_size, pressure_size>(pressure_index,
1414 .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
1415 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1416
1417 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
1418 dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
1419
1420 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1421 {
1422 double const alpha_bar =
1423 this->process_data_.micro_porosity_parameters
1424 ->mass_exchange_coefficient;
1425 auto const p_L_m =
1426 *std::get<MicroPressure>(this->current_states_[ip]);
1427 local_rhs.template segment<pressure_size>(pressure_index)
1428 .noalias() -=
1429 N_p.transpose() * alpha_bar / mu * (-p_cap_ip - p_L_m) * w;
1430
1431 local_Jac
1432 .template block<pressure_size, pressure_size>(pressure_index,
1434 .noalias() += N_p.transpose() * alpha_bar / mu * N_p * w;
1435 if (p_cap_ip != p_cap_prev_ip)
1436 {
1437 auto const p_L_m_prev = **std::get<PrevState<MicroPressure>>(
1438 this->prev_states_[ip]);
1439 local_Jac
1440 .template block<pressure_size, pressure_size>(
1442 .noalias() += N_p.transpose() * alpha_bar / mu *
1443 (p_L_m - p_L_m_prev) /
1444 (p_cap_ip - p_cap_prev_ip) * N_p * w;
1445 }
1446 }
1447 }
1448
1449 if (this->process_data_.apply_mass_lumping)
1450 {
1451 storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
1452 storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
1453 storage_p_a_S_Jpp =
1454 storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
1455 }
1456
1457 // pressure equation, pressure part.
1458 local_Jac
1459 .template block<pressure_size, pressure_size>(pressure_index,
1461 .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
1462
1463 // pressure equation, displacement part.
1464 local_Jac
1465 .template block<pressure_size, displacement_size>(pressure_index,
1467 .noalias() = Kpu / dt;
1468
1469 // pressure equation
1470 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
1471 laplace_p * p_L +
1472 (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
1473 Kpu * (u - u_prev) / dt;
1474
1475 // displacement equation
1476 local_rhs.template segment<displacement_size>(displacement_index)
1477 .noalias() += Kup * p_L;
1478}
static void assembleWithJacobianEvalConstitutiveSetting(double const t, double const dt, ParameterLib::SpatialPosition const &x_position, IpData &ip_data, MPL::VariableArray &variables, MPL::VariableArray &variables_prev, MPL::Medium const *const medium, TemperatureData const T_data, CapillaryPressureData< DisplacementDim > const &p_cap_data, ConstitutiveData< DisplacementDim > &CD, StatefulData< DisplacementDim > &SD, StatefulDataPrev< DisplacementDim > const &SD_prev, std::optional< MicroPorosityParameters > const &micro_porosity_parameters, MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material, ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > &material_state_data)
BaseLib::StrongType< double, struct TemperatureDataTag > TemperatureData
Definition Base.h:67
ConstitutiveModels< DisplacementDim > createConstitutiveModels(TRMProcessData const &process_data, MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
std::tuple< StiffnessTensor< DisplacementDim >, ProcessLib::ThermoRichardsMechanics::PorosityData, Density, LiquidDensity, ProcessLib::ThermoRichardsMechanics::BiotData, ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv, ProcessLib::ThermoRichardsMechanics::LiquidViscosityData, ProcessLib::ThermoRichardsMechanics::SolidCompressibilityData, ProcessLib::ThermoRichardsMechanics::BishopsData, PrevState< ProcessLib::ThermoRichardsMechanics::BishopsData >, ProcessLib::ThermoRichardsMechanics::PermeabilityData< DisplacementDim >, SaturationSecantDerivative > ConstitutiveData
Data that is needed for the equation system assembly.
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
VectorType< GlobalDim > GlobalDimVectorType

References MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), ProcessLib::RichardsMechanics::createConstitutiveModels(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::gas_phase_pressure, NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::temperature, and MathLib::KelvinVector::tensorToKelvin().

◆ assembleWithJacobianEvalConstitutiveSetting()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianEvalConstitutiveSetting ( double const t,
double const dt,
ParameterLib::SpatialPosition const & x_position,
IpData & ip_data,
MPL::VariableArray & variables,
MPL::VariableArray & variables_prev,
MPL::Medium const *const medium,
TemperatureData const T_data,
CapillaryPressureData< DisplacementDim > const & p_cap_data,
ConstitutiveData< DisplacementDim > & CD,
StatefulData< DisplacementDim > & SD,
StatefulDataPrev< DisplacementDim > const & SD_prev,
std::optional< MicroPorosityParameters > const & micro_porosity_parameters,
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material,
ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > & material_state_data )
staticprivate

Definition at line 754 of file RichardsMechanicsFEM-impl.h.

772{
773 auto const& liquid_phase = medium->phase("AqueousLiquid");
774 auto const& solid_phase = medium->phase("Solid");
775
776 auto const& identity2 = MathLib::KelvinVector::Invariants<
778 DisplacementDim)>::identity2;
779
780 double const temperature = T_data();
781 double const p_cap_ip = p_cap_data.p_cap;
782 double const p_cap_prev_ip = p_cap_data.p_cap_prev;
783
784 auto const& eps = std::get<StrainData<DisplacementDim>>(SD);
785 auto& S_L =
786 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(SD).S_L;
787 auto const S_L_prev =
788 std::get<
789 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
790 SD_prev)
791 ->S_L;
792 auto const alpha =
793 medium->property(MPL::PropertyType::biot_coefficient)
794 .template value<double>(variables, x_position, t, dt);
795 *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD) = alpha;
796
797 auto const C_el = ip_data.computeElasticTangentStiffness(
798 t, x_position, dt, temperature, solid_material,
799 *material_state_data.material_state_variables);
800
801 auto const beta_SR =
802 (1 - alpha) / solid_material.getBulkModulus(t, x_position, &C_el);
803 variables.grain_compressibility = beta_SR;
804 std::get<ProcessLib::ThermoRichardsMechanics::SolidCompressibilityData>(CD)
805 .beta_SR = beta_SR;
806
807 auto const rho_LR =
808 liquid_phase.property(MPL::PropertyType::density)
809 .template value<double>(variables, x_position, t, dt);
810 variables.density = rho_LR;
811 *std::get<LiquidDensity>(CD) = rho_LR;
812
813 S_L = medium->property(MPL::PropertyType::saturation)
814 .template value<double>(variables, x_position, t, dt);
815 variables.liquid_saturation = S_L;
816 variables_prev.liquid_saturation = S_L_prev;
817
818 // tangent derivative for Jacobian
819 double const dS_L_dp_cap =
820 medium->property(MPL::PropertyType::saturation)
821 .template dValue<double>(variables,
822 MPL::Variable::capillary_pressure,
823 x_position, t, dt);
824 std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(CD)
825 .dS_L_dp_cap = dS_L_dp_cap;
826 // secant derivative from time discretization for storage
827 // use tangent, if secant is not available
828 double const DeltaS_L_Deltap_cap =
829 (p_cap_ip == p_cap_prev_ip)
830 ? dS_L_dp_cap
831 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
832 std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap =
833 DeltaS_L_Deltap_cap;
834
835 auto const chi = [medium, x_position, t, dt](double const S_L)
836 {
838 vs.liquid_saturation = S_L;
839 return medium->property(MPL::PropertyType::bishops_effective_stress)
840 .template value<double>(vs, x_position, t, dt);
841 };
842 double const chi_S_L = chi(S_L);
843 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD).chi_S_L =
844 chi_S_L;
845 double const chi_S_L_prev = chi(S_L_prev);
846 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::BishopsData>>(CD)
847 ->chi_S_L = chi_S_L_prev;
848
849 auto const dchi_dS_L =
850 medium->property(MPL::PropertyType::bishops_effective_stress)
851 .template dValue<double>(
852 variables, MPL::Variable::liquid_saturation, x_position, t, dt);
853 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD).dchi_dS_L =
854 dchi_dS_L;
855
856 double const p_FR = -chi_S_L * p_cap_ip;
857 variables.effective_pore_pressure = p_FR;
858 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
859
860 // Set volumetric strain rate for the general case without swelling.
861 variables.volumetric_strain = Invariants::trace(eps.eps);
862 // TODO (CL) changed that, using eps_prev for the moment, not B * u_prev
863 // variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
864 variables_prev.volumetric_strain = Invariants::trace(
865 std::get<PrevState<StrainData<DisplacementDim>>>(SD_prev)->eps);
866
867 auto& phi =
868 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(SD).phi;
869 { // Porosity update
870 auto const phi_prev =
871 std::get<
872 PrevState<ProcessLib::ThermoRichardsMechanics::PorosityData>>(
873 SD_prev)
874 ->phi;
875 variables_prev.porosity = phi_prev;
876 phi = medium->property(MPL::PropertyType::porosity)
877 .template value<double>(variables, variables_prev, x_position,
878 t, dt);
879 variables.porosity = phi;
880 }
881 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi = phi;
882
883 if (alpha < phi)
884 {
885 auto const eid =
886 x_position.getElementID()
887 ? static_cast<std::ptrdiff_t>(*x_position.getElementID())
888 : static_cast<std::ptrdiff_t>(-1);
889 auto const ip =
890 x_position.getIntegrationPoint()
891 ? static_cast<std::ptrdiff_t>(*x_position.getIntegrationPoint())
892 : static_cast<std::ptrdiff_t>(-1);
893 OGS_FATAL(
894 "RichardsMechanics: Biot-coefficient {} is smaller than "
895 "porosity {} in element/integration point {}/{}.",
896 alpha, phi, eid, ip);
897 }
898
899 auto const mu = liquid_phase.property(MPL::PropertyType::viscosity)
900 .template value<double>(variables, x_position, t, dt);
901 *std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(CD) =
902 mu;
903
904 {
905 // Swelling and possibly volumetric strain rate update.
906 auto& sigma_sw =
907 std::get<ProcessLib::ThermoRichardsMechanics::
908 ConstitutiveStress_StrainTemperature::
909 SwellingDataStateful<DisplacementDim>>(SD);
910 auto const& sigma_sw_prev =
911 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
912 ConstitutiveStress_StrainTemperature::
913 SwellingDataStateful<DisplacementDim>>>(
914 SD_prev);
915 auto const transport_porosity_prev = std::get<PrevState<
917 SD_prev);
918 auto const phi_prev = std::get<
919 PrevState<ProcessLib::ThermoRichardsMechanics::PorosityData>>(
920 SD_prev);
921 auto& transport_porosity = std::get<
923 auto& p_L_m = std::get<MicroPressure>(SD);
924 auto const p_L_m_prev = std::get<PrevState<MicroPressure>>(SD_prev);
925 auto& S_L_m = std::get<MicroSaturation>(SD);
926 auto const S_L_m_prev = std::get<PrevState<MicroSaturation>>(SD_prev);
927
929 *medium, solid_phase, C_el, rho_LR, mu, micro_porosity_parameters,
930 alpha, phi, p_cap_ip, variables, variables_prev, x_position, t, dt,
931 sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev,
932 transport_porosity, p_L_m_prev, S_L_m_prev, p_L_m, S_L_m);
933 }
934
935 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
936 {
937 if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
938 {
939 auto& transport_porosity =
940 std::get<
942 SD)
943 .phi;
944 auto const transport_porosity_prev = std::get<PrevState<
946 SD_prev)
947 ->phi;
948 variables_prev.transport_porosity = transport_porosity_prev;
949
951 medium->property(MPL::PropertyType::transport_porosity)
952 .template value<double>(variables, variables_prev,
953 x_position, t, dt);
955 }
956 }
957 else
958 {
959 variables.transport_porosity = phi;
960 }
961
962 // Set mechanical variables for the intrinsic permeability model
963 // For stress dependent permeability.
964 {
965 // TODO mechanical constitutive relation will be evaluated afterwards
966 auto const sigma_total =
967 (std::get<ProcessLib::ThermoRichardsMechanics::
968 ConstitutiveStress_StrainTemperature::
969 EffectiveStressData<DisplacementDim>>(SD)
970 .sigma_eff +
971 alpha * p_FR * identity2)
972 .eval();
973 // For stress dependent permeability.
974 variables.total_stress.emplace<SymmetricTensor>(
976 }
977
978 variables.equivalent_plastic_strain =
979 material_state_data.material_state_variables
980 ->getEquivalentPlasticStrain();
981
982 double const k_rel =
983 medium->property(MPL::PropertyType::relative_permeability)
984 .template value<double>(variables, x_position, t, dt);
985
986 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
987 medium->property(MPL::PropertyType::permeability)
988 .value(variables, x_position, t, dt));
989
990 std::get<
992 CD)
993 .k_rel = k_rel;
994 std::get<
996 CD)
997 .Ki = K_intrinsic;
998
999 //
1000 // displacement equation, displacement part
1001 //
1002
1003 {
1004 auto& sigma_sw =
1005 std::get<ProcessLib::ThermoRichardsMechanics::
1006 ConstitutiveStress_StrainTemperature::
1007 SwellingDataStateful<DisplacementDim>>(SD)
1008 .sigma_sw;
1009
1010 auto& eps_m =
1011 std::get<ProcessLib::ThermoRichardsMechanics::
1012 ConstitutiveStress_StrainTemperature::
1013 MechanicalStrainData<DisplacementDim>>(SD)
1014 .eps_m;
1015 eps_m.noalias() =
1016 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
1017 ? eps.eps + C_el.inverse() * sigma_sw
1018 : eps.eps;
1019 variables.mechanical_strain
1021 eps_m);
1022 }
1023
1024 {
1025 auto& sigma_eff =
1026 std::get<ProcessLib::ThermoRichardsMechanics::
1027 ConstitutiveStress_StrainTemperature::
1028 EffectiveStressData<DisplacementDim>>(SD);
1029 auto const& sigma_eff_prev =
1030 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
1031 ConstitutiveStress_StrainTemperature::
1032 EffectiveStressData<DisplacementDim>>>(
1033 SD_prev);
1034 auto const& eps_m =
1035 std::get<ProcessLib::ThermoRichardsMechanics::
1036 ConstitutiveStress_StrainTemperature::
1037 MechanicalStrainData<DisplacementDim>>(SD);
1038 auto& eps_m_prev =
1039 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
1040 ConstitutiveStress_StrainTemperature::
1041 MechanicalStrainData<DisplacementDim>>>(
1042 SD_prev);
1043
1044 auto C = ip_data.updateConstitutiveRelation(
1045 variables, t, x_position, dt, temperature, sigma_eff,
1046 sigma_eff_prev, eps_m, eps_m_prev, solid_material,
1047 material_state_data.material_state_variables);
1048
1049 *std::get<StiffnessTensor<DisplacementDim>>(CD) = std::move(C);
1050 }
1051
1052 // p_SR
1053 variables.solid_grain_pressure =
1054 p_FR -
1055 std::get<ProcessLib::ThermoRichardsMechanics::
1056 ConstitutiveStress_StrainTemperature::EffectiveStressData<
1057 DisplacementDim>>(SD)
1058 .sigma_eff.dot(identity2) /
1059 (3 * (1 - phi));
1060 auto const rho_SR =
1061 solid_phase.property(MPL::PropertyType::density)
1062 .template value<double>(variables, x_position, t, dt);
1063
1064 double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
1065 *std::get<Density>(CD) = rho;
1066}
void updateSwellingStressAndVolumetricStrain(MaterialPropertyLib::Medium const &medium, MaterialPropertyLib::Phase const &solid_phase, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > const &C_el, double const rho_LR, double const mu, std::optional< MicroPorosityParameters > micro_porosity_parameters, double const alpha, double const phi, double const p_cap_ip, MPL::VariableArray &variables, MPL::VariableArray &variables_prev, ParameterLib::SpatialPosition const &x_position, double const t, double const dt, ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature::SwellingDataStateful< DisplacementDim > &sigma_sw, PrevState< ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature::SwellingDataStateful< DisplacementDim > > const &sigma_sw_prev, PrevState< ProcessLib::ThermoRichardsMechanics::TransportPorosityData > const phi_M_prev, PrevState< ProcessLib::ThermoRichardsMechanics::PorosityData > const phi_prev, ProcessLib::ThermoRichardsMechanics::TransportPorosityData &phi_M, PrevState< MicroPressure > const p_L_m_prev, PrevState< MicroSaturation > const S_L_m_prev, MicroPressure &p_L_m, MicroSaturation &S_L_m)
virtual double getBulkModulus(double const, ParameterLib::SpatialPosition const &, KelvinMatrix const *const =nullptr) const

References MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialLib::Solids::MechanicsBase< DisplacementDim >::getBulkModulus(), ParameterLib::SpatialPosition::getElementID(), ParameterLib::SpatialPosition::getIntegrationPoint(), MaterialPropertyLib::VariableArray::grain_compressibility, MaterialPropertyLib::Medium::hasProperty(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_saturation, ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim >::material_state_variables, MaterialPropertyLib::VariableArray::mechanical_strain, OGS_FATAL, ProcessLib::RichardsMechanics::CapillaryPressureData< DisplacementDim >::p_cap, ProcessLib::RichardsMechanics::CapillaryPressureData< DisplacementDim >::p_cap_prev, MaterialPropertyLib::Medium::phase(), MaterialPropertyLib::VariableArray::porosity, MaterialPropertyLib::Medium::property(), MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::VariableArray::transport_porosity, ProcessLib::RichardsMechanics::updateSwellingStressAndVolumetricStrain(), MaterialPropertyLib::Property::value(), and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ assembleWithJacobianForDeformationEquations()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForDeformationEquations ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev,
std::vector< double > & local_b_data,
std::vector< double > & local_Jac_data )
private

Assemble local matrices and vectors arise from the linearized discretized weak form of the residual of the momentum balance equation,

\nabla (\sigma - \alpha_b p \mathrm{I}) = f

where \sigma is the effective stress tensor, p is the pore pressure, \alpha_b is the Biot constant, \mathrm{I} is the identity tensor, and f is the body force.

Parameters
tTime
dtTime increment
local_xNodal values of x of an element.
local_x_prevNodal values of x_{prev} of an element.
local_b_dataRight hand side vector of an element.
local_Jac_dataElement Jacobian matrix for the Newton-Raphson method.

Definition at line 1497 of file RichardsMechanicsFEM-impl.h.

1504{
1505 OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
1506}

References OGS_FATAL.

◆ assembleWithJacobianForPressureEquations()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForPressureEquations ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev,
std::vector< double > & local_b_data,
std::vector< double > & local_Jac_data )
private

Assemble local matrices and vectors arise from the linearized discretized weak form of the residual of the mass balance equation of single phase flow,

\alpha \cdot{p} - \nabla (K (\nabla p + \rho g \nabla z) + \alpha_b \nabla \cdot \dot{u} = Q

where alpha is a coefficient may depend on storage or the fluid density change, \rho is the fluid density, g is the gravitational acceleration, z is the vertical coordinate, u is the displacement, and Q is the source/sink term.

Parameters
tTime
dtTime increment
local_xNodal values of x of an element.
local_x_prevNodal values of x_{prev} of an element.
local_b_dataRight hand side vector of an element.
local_Jac_dataElement Jacobian matrix for the Newton-Raphson method.

Definition at line 1483 of file RichardsMechanicsFEM-impl.h.

1490{
1491 OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
1492}

References OGS_FATAL.

◆ assembleWithJacobianForStaggeredScheme()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::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 )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1511 of file RichardsMechanicsFEM-impl.h.

1518{
1519 // For the equations with pressure
1520 if (process_id == 0)
1521 {
1522 assembleWithJacobianForPressureEquations(t, dt, local_x, local_x_prev,
1523 local_b_data, local_Jac_data);
1524 return;
1525 }
1526
1527 // For the equations with deformation
1528 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_x_prev,
1529 local_b_data, local_Jac_data);
1530}
void assembleWithJacobianForPressureEquations(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void assembleWithJacobianForDeformationEquations(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)

◆ computeSecondaryVariableConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< 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 1535 of file RichardsMechanicsFEM-impl.h.

1539{
1540 auto p_L = local_x.template segment<pressure_size>(pressure_index);
1541 auto u = local_x.template segment<displacement_size>(displacement_index);
1542
1543 auto p_L_prev =
1544 local_x_prev.template segment<pressure_size>(pressure_index);
1545 auto u_prev =
1546 local_x_prev.template segment<displacement_size>(displacement_index);
1547
1548 auto const& identity2 = MathLib::KelvinVector::Invariants<
1550 DisplacementDim)>::identity2;
1551
1552 auto const& medium =
1553 this->process_data_.media_map.getMedium(this->element_.getID());
1554 auto const& liquid_phase = medium->phase("AqueousLiquid");
1555 auto const& solid_phase = medium->phase("Solid");
1556 MPL::VariableArray variables;
1557 MPL::VariableArray variables_prev;
1558
1560 x_position.setElementID(this->element_.getID());
1561
1562 unsigned const n_integration_points =
1564
1565 double saturation_avg = 0;
1566 double porosity_avg = 0;
1567
1569 KV sigma_avg = KV::Zero();
1570
1571 for (unsigned ip = 0; ip < n_integration_points; ip++)
1572 {
1573 x_position.setIntegrationPoint(ip);
1574 auto const& N_p = _ip_data[ip].N_p;
1575 auto const& N_u = _ip_data[ip].N_u;
1576 auto const& dNdx_u = _ip_data[ip].dNdx_u;
1577
1578 auto const x_coord =
1579 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1581 this->element_, N_u);
1582 auto const B =
1583 LinearBMatrix::computeBMatrix<DisplacementDim,
1584 ShapeFunctionDisplacement::NPOINTS,
1586 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
1587
1588 double p_cap_ip;
1589 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
1590
1591 double p_cap_prev_ip;
1592 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
1593
1594 variables.capillary_pressure = p_cap_ip;
1595 variables.liquid_phase_pressure = -p_cap_ip;
1596 // setting pG to 1 atm
1597 // TODO : rewrite equations s.t. p_L = pG-p_cap
1598 variables.gas_phase_pressure = 1.0e5;
1599
1600 auto const temperature =
1601 medium->property(MPL::PropertyType::reference_temperature)
1602 .template value<double>(variables, x_position, t, dt);
1603 variables.temperature = temperature;
1604
1605 auto& eps =
1606 std::get<StrainData<DisplacementDim>>(this->current_states_[ip])
1607 .eps;
1608 eps.noalias() = B * u;
1609 auto& S_L =
1610 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
1611 this->current_states_[ip])
1612 .S_L;
1613 auto const S_L_prev =
1614 std::get<
1615 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
1616 this->prev_states_[ip])
1617 ->S_L;
1618 S_L = medium->property(MPL::PropertyType::saturation)
1619 .template value<double>(variables, x_position, t, dt);
1620 variables.liquid_saturation = S_L;
1621 variables_prev.liquid_saturation = S_L_prev;
1622
1623 auto const chi = [medium, x_position, t, dt](double const S_L)
1624 {
1626 vs.liquid_saturation = S_L;
1627 return medium->property(MPL::PropertyType::bishops_effective_stress)
1628 .template value<double>(vs, x_position, t, dt);
1629 };
1630 double const chi_S_L = chi(S_L);
1631 double const chi_S_L_prev = chi(S_L_prev);
1632
1633 auto const alpha =
1634 medium->property(MPL::PropertyType::biot_coefficient)
1635 .template value<double>(variables, x_position, t, dt);
1636
1637 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
1638 t, x_position, dt, temperature, this->solid_material_,
1639 *this->material_states_[ip].material_state_variables);
1640
1641 auto const beta_SR = (1 - alpha) / this->solid_material_.getBulkModulus(
1642 t, x_position, &C_el);
1643 variables.grain_compressibility = beta_SR;
1644
1645 variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
1646 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
1647
1648 // Set volumetric strain rate for the general case without swelling.
1649 variables.volumetric_strain = Invariants::trace(eps);
1650 variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
1651
1652 auto& phi = std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
1653 this->current_states_[ip])
1654 .phi;
1655 { // Porosity update
1656 auto const phi_prev = std::get<PrevState<
1658 this->prev_states_[ip])
1659 ->phi;
1660 variables_prev.porosity = phi_prev;
1661 phi = medium->property(MPL::PropertyType::porosity)
1662 .template value<double>(variables, variables_prev,
1663 x_position, t, dt);
1664 variables.porosity = phi;
1665 }
1666
1667 auto const rho_LR =
1668 liquid_phase.property(MPL::PropertyType::density)
1669 .template value<double>(variables, x_position, t, dt);
1670 variables.density = rho_LR;
1671 auto const mu =
1672 liquid_phase.property(MPL::PropertyType::viscosity)
1673 .template value<double>(variables, x_position, t, dt);
1674
1675 {
1676 // Swelling and possibly volumetric strain rate update.
1677 auto& sigma_sw =
1678 std::get<ProcessLib::ThermoRichardsMechanics::
1679 ConstitutiveStress_StrainTemperature::
1680 SwellingDataStateful<DisplacementDim>>(
1681 this->current_states_[ip]);
1682 auto const& sigma_sw_prev = std::get<
1683 PrevState<ProcessLib::ThermoRichardsMechanics::
1684 ConstitutiveStress_StrainTemperature::
1685 SwellingDataStateful<DisplacementDim>>>(
1686 this->prev_states_[ip]);
1687 auto const transport_porosity_prev = std::get<PrevState<
1689 this->prev_states_[ip]);
1690 auto const phi_prev = std::get<
1691 PrevState<ProcessLib::ThermoRichardsMechanics::PorosityData>>(
1692 this->prev_states_[ip]);
1693 auto& transport_porosity = std::get<
1695 this->current_states_[ip]);
1696 auto& p_L_m = std::get<MicroPressure>(this->current_states_[ip]);
1697 auto const p_L_m_prev =
1698 std::get<PrevState<MicroPressure>>(this->prev_states_[ip]);
1699 auto& S_L_m = std::get<MicroSaturation>(this->current_states_[ip]);
1700 auto const S_L_m_prev =
1701 std::get<PrevState<MicroSaturation>>(this->prev_states_[ip]);
1702
1704 *medium, solid_phase, C_el, rho_LR, mu,
1705 this->process_data_.micro_porosity_parameters, alpha, phi,
1706 p_cap_ip, variables, variables_prev, x_position, t, dt,
1707 sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev,
1708 transport_porosity, p_L_m_prev, S_L_m_prev, p_L_m, S_L_m);
1709 }
1710
1711 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
1712 {
1713 if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
1714 {
1715 auto& transport_porosity =
1716 std::get<ProcessLib::ThermoRichardsMechanics::
1717 TransportPorosityData>(
1718 this->current_states_[ip])
1719 .phi;
1720 auto const transport_porosity_prev =
1721 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
1722 TransportPorosityData>>(
1723 this->prev_states_[ip])
1724 ->phi;
1725
1726 variables_prev.transport_porosity = transport_porosity_prev;
1727
1729 medium->property(MPL::PropertyType::transport_porosity)
1730 .template value<double>(variables, variables_prev,
1731 x_position, t, dt);
1733 }
1734 }
1735 else
1736 {
1737 variables.transport_porosity = phi;
1738 }
1739
1740 auto const& sigma_eff =
1741 std::get<ProcessLib::ThermoRichardsMechanics::
1742 ConstitutiveStress_StrainTemperature::
1743 EffectiveStressData<DisplacementDim>>(
1744 this->current_states_[ip])
1745 .sigma_eff;
1746
1747 // Set mechanical variables for the intrinsic permeability model
1748 // For stress dependent permeability.
1749 {
1750 auto const sigma_total =
1751 (sigma_eff + alpha * chi_S_L * identity2 * p_cap_ip).eval();
1752 // For stress dependent permeability.
1753 variables.total_stress.emplace<SymmetricTensor>(
1755 sigma_total));
1756 }
1757
1758 variables.equivalent_plastic_strain =
1759 this->material_states_[ip]
1760 .material_state_variables->getEquivalentPlasticStrain();
1761
1762 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
1763 medium->property(MPL::PropertyType::permeability)
1764 .value(variables, x_position, t, dt));
1765
1766 double const k_rel =
1767 medium->property(MPL::PropertyType::relative_permeability)
1768 .template value<double>(variables, x_position, t, dt);
1769
1770 GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
1771
1772 double const p_FR = -chi_S_L * p_cap_ip;
1773 // p_SR
1774 variables.solid_grain_pressure =
1775 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
1776 auto const rho_SR =
1777 solid_phase.property(MPL::PropertyType::density)
1778 .template value<double>(variables, x_position, t, dt);
1779 *std::get<DrySolidDensity>(this->output_data_[ip]) = (1 - phi) * rho_SR;
1780
1781 {
1782 auto& SD = this->current_states_[ip];
1783 auto const& sigma_sw =
1784 std::get<ProcessLib::ThermoRichardsMechanics::
1785 ConstitutiveStress_StrainTemperature::
1786 SwellingDataStateful<DisplacementDim>>(SD)
1787 .sigma_sw;
1788 auto& eps_m =
1789 std::get<ProcessLib::ThermoRichardsMechanics::
1790 ConstitutiveStress_StrainTemperature::
1791 MechanicalStrainData<DisplacementDim>>(SD)
1792 .eps_m;
1793 eps_m.noalias() =
1794 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
1795 ? eps + C_el.inverse() * sigma_sw
1796 : eps;
1797 variables.mechanical_strain.emplace<
1799 eps_m);
1800 }
1801
1802 {
1803 auto& SD = this->current_states_[ip];
1804 auto const& SD_prev = this->prev_states_[ip];
1805 auto& sigma_eff =
1806 std::get<ProcessLib::ThermoRichardsMechanics::
1807 ConstitutiveStress_StrainTemperature::
1808 EffectiveStressData<DisplacementDim>>(SD);
1809 auto const& sigma_eff_prev = std::get<
1810 PrevState<ProcessLib::ThermoRichardsMechanics::
1811 ConstitutiveStress_StrainTemperature::
1812 EffectiveStressData<DisplacementDim>>>(
1813 SD_prev);
1814 auto const& eps_m =
1815 std::get<ProcessLib::ThermoRichardsMechanics::
1816 ConstitutiveStress_StrainTemperature::
1817 MechanicalStrainData<DisplacementDim>>(SD);
1818 auto const& eps_m_prev = std::get<
1819 PrevState<ProcessLib::ThermoRichardsMechanics::
1820 ConstitutiveStress_StrainTemperature::
1821 MechanicalStrainData<DisplacementDim>>>(
1822 SD_prev);
1823
1824 _ip_data[ip].updateConstitutiveRelation(
1825 variables, t, x_position, dt, temperature, sigma_eff,
1826 sigma_eff_prev, eps_m, eps_m_prev, this->solid_material_,
1827 this->material_states_[ip].material_state_variables);
1828 }
1829
1830 auto const& b = this->process_data_.specific_body_force;
1831
1832 // Compute the velocity
1833 auto const& dNdx_p = _ip_data[ip].dNdx_p;
1834 std::get<
1836 this->output_data_[ip])
1837 ->noalias() = -K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b;
1838
1839 saturation_avg += S_L;
1840 porosity_avg += phi;
1841 sigma_avg += sigma_eff;
1842 }
1843 saturation_avg /= n_integration_points;
1844 porosity_avg /= n_integration_points;
1845 sigma_avg /= n_integration_points;
1846
1847 (*this->process_data_.element_saturation)[this->element_.getID()] =
1848 saturation_avg;
1849 (*this->process_data_.element_porosity)[this->element_.getID()] =
1850 porosity_avg;
1851
1852 Eigen::Map<KV>(
1853 &(*this->process_data_.element_stresses)[this->element_.getID() *
1854 KV::RowsAtCompileTime]) =
1856
1858 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1859 DisplacementDim>(this->element_, this->is_axially_symmetric_, p_L,
1860 *this->process_data_.pressure_interpolated);
1861}
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
std::vector< OutputData< DisplacementDim > > output_data_

References MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::VariableArray::grain_compressibility, NumLib::interpolateToHigherOrderNodes(), NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::VariableArray::porosity, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::VariableArray::transport_porosity, ProcessLib::RichardsMechanics::updateSwellingStressAndVolumetricStrain(), and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ getShapeMatrix()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
Eigen::Map< const Eigen::RowVectorXd > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< 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 165 of file RichardsMechanicsFEM.h.

167 {
168 auto const& N_u = _secondary_data.N_u[integration_point];
169
170 // assumes N is stored contiguously in memory
171 return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
172 }

References ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_secondary_data, and ProcessLib::RichardsMechanics::SecondaryData< ShapeMatrixType >::N_u.

◆ initializeConcrete()

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

Set initial stress from parameter.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 116 of file RichardsMechanicsFEM.h.

117 {
118 unsigned const n_integration_points =
120
121 for (unsigned ip = 0; ip < n_integration_points; ip++)
122 {
123 auto& SD = this->current_states_[ip];
124 auto& ip_data = _ip_data[ip];
125
126 ParameterLib::SpatialPosition const x_position{
127 std::nullopt, this->element_.getID(), ip,
129 ShapeFunctionDisplacement,
131 ip_data.N_u))};
132
134 if (this->process_data_.initial_stress.value)
135 {
136 std::get<ProcessLib::ThermoRichardsMechanics::
137 ConstitutiveStress_StrainTemperature::
138 EffectiveStressData<DisplacementDim>>(SD)
139 .sigma_eff =
141 DisplacementDim>(
142 // The data in process_data_.initial_stress.value can
143 // be total stress or effective stress.
144 (*this->process_data_.initial_stress.value)(
145 std::numeric_limits<
146 double>::quiet_NaN() /* time independent */,
147 x_position));
148 }
149
150 double const t = 0; // TODO (naumov) pass t from top
151 this->solid_material_.initializeInternalStateVariables(
152 t, x_position,
153 *this->material_states_[ip].material_state_variables);
154
155 this->material_states_[ip].pushBackState();
156
157 this->prev_states_[ip] = SD;
158 }
159 }
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)

References ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::current_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::element_, MeshLib::Element::getID(), NumLib::GenericIntegrationMethod::getNumberOfPoints(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::material_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::prev_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::process_data_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::solid_material_, and MathLib::KelvinVector::symmetricTensorToKelvinVector().

◆ setInitialConditionsConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 203 of file RichardsMechanicsFEM-impl.h.

207{
208 assert(local_x.size() == pressure_size + displacement_size);
209
210 auto const p_L = local_x.template segment<pressure_size>(pressure_index);
211
212 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
213 auto const& medium =
214 this->process_data_.media_map.getMedium(this->element_.getID());
215 MPL::VariableArray variables;
216
218 x_position.setElementID(this->element_.getID());
219
220 auto const& solid_phase = medium->phase("Solid");
221
222 auto const& identity2 = MathLib::KelvinVector::Invariants<
224 DisplacementDim)>::identity2;
225
226 unsigned const n_integration_points =
228 for (unsigned ip = 0; ip < n_integration_points; ip++)
229 {
230 x_position.setIntegrationPoint(ip);
231
232 auto const& N_p = _ip_data[ip].N_p;
233
234 double p_cap_ip;
235 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
236
237 variables.capillary_pressure = p_cap_ip;
238 variables.liquid_phase_pressure = -p_cap_ip;
239 // setting pG to 1 atm
240 // TODO : rewrite equations s.t. p_L = pG-p_cap
241 variables.gas_phase_pressure = 1.0e5;
242
243 {
244 auto& p_L_m = std::get<MicroPressure>(this->current_states_[ip]);
245 auto& p_L_m_prev =
246 std::get<PrevState<MicroPressure>>(this->prev_states_[ip]);
247 **p_L_m_prev = -p_cap_ip;
248 *p_L_m = -p_cap_ip;
249 }
250
251 auto const temperature =
252 medium->property(MPL::PropertyType::reference_temperature)
253 .template value<double>(variables, x_position, t, dt);
254 variables.temperature = temperature;
255
256 auto& S_L_prev =
257 std::get<
258 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
259 this->prev_states_[ip])
260 ->S_L;
261 S_L_prev = medium->property(MPL::PropertyType::saturation)
262 .template value<double>(variables, x_position, t, dt);
263
264 if (this->process_data_.initial_stress.isTotalStress())
265 {
266 auto const alpha_b =
267 medium->property(MPL::PropertyType::biot_coefficient)
268 .template value<double>(variables, x_position, t, dt);
269
270 variables.liquid_saturation = S_L_prev;
271 double const chi_S_L =
272 medium->property(MPL::PropertyType::bishops_effective_stress)
273 .template value<double>(variables, x_position, t, dt);
274
275 // Initial stresses are total stress, which were assigned to
276 // sigma_eff in
277 // RichardsMechanicsLocalAssembler::initializeConcrete().
278 auto& sigma_eff =
279 std::get<ProcessLib::ThermoRichardsMechanics::
280 ConstitutiveStress_StrainTemperature::
281 EffectiveStressData<DisplacementDim>>(
282 this->current_states_[ip]);
283
284 auto& sigma_eff_prev = std::get<
285 PrevState<ProcessLib::ThermoRichardsMechanics::
286 ConstitutiveStress_StrainTemperature::
287 EffectiveStressData<DisplacementDim>>>(
288 this->prev_states_[ip]);
289
290 // Reset sigma_eff to effective stress
291 sigma_eff.sigma_eff.noalias() +=
292 chi_S_L * alpha_b * (-p_cap_ip) * identity2;
293 sigma_eff_prev->sigma_eff = sigma_eff.sigma_eff;
294 }
295
296 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
297 {
299 vars.capillary_pressure = p_cap_ip;
300
301 auto& S_L_m = std::get<MicroSaturation>(this->current_states_[ip]);
302 auto& S_L_m_prev =
303 std::get<PrevState<MicroSaturation>>(this->prev_states_[ip]);
304
305 *S_L_m = medium->property(MPL::PropertyType::saturation_micro)
306 .template value<double>(vars, x_position, t, dt);
307 *S_L_m_prev = S_L_m;
308 }
309
310 // Set eps_m_prev from potentially non-zero eps and sigma_sw from
311 // restart.
312 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
313 t, x_position, dt, temperature, this->solid_material_,
314 *this->material_states_[ip].material_state_variables);
315 auto& eps =
316 std::get<StrainData<DisplacementDim>>(this->current_states_[ip])
317 .eps;
318 auto const& sigma_sw =
319 std::get<ProcessLib::ThermoRichardsMechanics::
320 ConstitutiveStress_StrainTemperature::
321 SwellingDataStateful<DisplacementDim>>(
322 this->current_states_[ip])
323 .sigma_sw;
324 auto& eps_m_prev =
325 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
326 ConstitutiveStress_StrainTemperature::
327 MechanicalStrainData<DisplacementDim>>>(
328 this->prev_states_[ip])
329 ->eps_m;
330
331 eps_m_prev.noalias() =
332 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
333 ? eps + C_el.inverse() * sigma_sw
334 : eps;
335 }
336}

References MaterialPropertyLib::VariableArray::capillary_pressure, MaterialPropertyLib::VariableArray::gas_phase_pressure, MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::detail::shapeFunctionInterpolate(), and MaterialPropertyLib::VariableArray::temperature.

Member Data Documentation

◆ _ip_data

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector<IpData, Eigen::aligned_allocator<IpData> > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data
private

◆ _secondary_data

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

◆ displacement_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::displacement_index = ShapeFunctionPressure::NPOINTS
staticprivate

Definition at line 248 of file RichardsMechanicsFEM.h.

◆ displacement_size

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

Definition at line 249 of file RichardsMechanicsFEM.h.

◆ KelvinVectorSize

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

Definition at line 72 of file RichardsMechanicsFEM.h.

◆ N_u_op

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
auto& ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< 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 78 of file RichardsMechanicsFEM.h.

◆ pressure_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
const int ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::pressure_index = 0
staticprivate

Definition at line 246 of file RichardsMechanicsFEM.h.

◆ pressure_size

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

Definition at line 247 of file RichardsMechanicsFEM.h.


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