129 ShapeFunctionPressure, DisplacementDim>::
130 RichardsMechanicsLocalAssembler(
134 bool const is_axially_symmetric,
137 e, integration_method, is_axially_symmetric, process_data}
139 unsigned const n_integration_points =
140 this->integration_method_.getNumberOfPoints();
142 ip_data_.resize(n_integration_points);
143 secondary_data_.N_u.resize(n_integration_points);
145 auto const shape_matrices_u =
147 ShapeMatricesTypeDisplacement,
148 DisplacementDim>(e, is_axially_symmetric,
149 this->integration_method_);
151 auto const shape_matrices_p =
153 ShapeMatricesTypePressure, DisplacementDim>(
154 e, is_axially_symmetric, this->integration_method_);
157 this->process_data_.media_map.getMedium(this->element_.getID());
159 for (
unsigned ip = 0; ip < n_integration_points; ip++)
161 auto& ip_data = ip_data_[ip];
162 auto const& sm_u = shape_matrices_u[ip];
163 ip_data_[ip].integration_weight =
164 this->integration_method_.getWeightedPoint(ip).getWeight() *
165 sm_u.integralMeasure * sm_u.detJ;
167 ip_data.N_u = sm_u.N;
168 ip_data.dNdx_u = sm_u.dNdx;
171 std::nullopt, this->element_.getID(),
174 ShapeMatricesTypeDisplacement>(
175 this->element_, ip_data.N_u))};
177 ip_data.N_p = shape_matrices_p[ip].N;
178 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
182 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
183 this->current_states_[ip])
186 .template initialValue<double>(
189 double>::quiet_NaN() );
191 auto& transport_porosity =
194 this->current_states_[ip])
196 transport_porosity = porosity;
197 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
201 .template initialValue<double>(
204 double>::quiet_NaN() );
207 secondary_data_.N_u[ip] = shape_matrices_u[ip].N;
214 ShapeFunctionPressure, DisplacementDim>::
215 setInitialConditionsConcrete(Eigen::VectorXd
const local_x,
219 assert(local_x.size() == pressure_size + displacement_size);
221 auto const [p_L, u] = localDOF(local_x);
223 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
225 this->process_data_.media_map.getMedium(this->element_.getID());
228 auto const& solid_phase = medium->phase(
"Solid");
232 DisplacementDim)>::identity2;
234 unsigned const n_integration_points =
235 this->integration_method_.getNumberOfPoints();
236 for (
unsigned ip = 0; ip < n_integration_points; ip++)
238 auto const& N_p = ip_data_[ip].N_p;
241 std::nullopt, this->element_.getID(),
245 this->element_, N_p))};
257 auto& p_L_m = std::get<MicroPressure>(this->current_states_[ip]);
259 std::get<PrevState<MicroPressure>>(this->prev_states_[ip]);
260 **p_L_m_prev = -p_cap_ip;
264 auto const temperature =
265 medium->property(MPL::PropertyType::reference_temperature)
266 .template value<double>(variables, x_position, t, dt);
272 this->prev_states_[ip])
274 S_L_prev = medium->property(MPL::PropertyType::saturation)
275 .template value<double>(variables, x_position, t, dt);
277 if (this->process_data_.initial_stress.isTotalStress())
280 medium->property(MPL::PropertyType::biot_coefficient)
281 .template value<double>(variables, x_position, t, dt);
284 double const chi_S_L =
285 medium->property(MPL::PropertyType::bishops_effective_stress)
286 .template value<double>(variables, x_position, t, dt);
293 DisplacementDim>>(this->current_states_[ip]);
295 auto& sigma_eff_prev =
296 std::get<
PrevState<ProcessLib::ConstitutiveRelations::
297 EffectiveStressData<DisplacementDim>>>(
298 this->prev_states_[ip]);
301 sigma_eff.sigma_eff.noalias() +=
302 chi_S_L * alpha_b * (-p_cap_ip) * identity2;
303 sigma_eff_prev->sigma_eff = sigma_eff.sigma_eff;
306 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
311 auto& S_L_m = std::get<MicroSaturation>(this->current_states_[ip]);
313 std::get<PrevState<MicroSaturation>>(this->prev_states_[ip]);
315 *S_L_m = medium->property(MPL::PropertyType::saturation_micro)
316 .template value<double>(vars, x_position, t, dt);
322 auto& SD = this->current_states_[ip];
325 DisplacementDim>>(SD)
328 auto const& N_u = ip_data_[ip].N_u;
329 auto const& dNdx_u = ip_data_[ip].dNdx_u;
333 this->element_, N_u);
336 ShapeFunctionDisplacement::NPOINTS,
338 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
340 std::get<StrainData<DisplacementDim>>(this->current_states_[ip])
342 eps.noalias() = B * u;
349 auto const C_el = ip_data_[ip].computeElasticTangentStiffness(
350 variables, t, x_position, dt, this->solid_material_,
351 *this->material_states_[ip].material_state_variables);
353 auto const& sigma_sw =
354 std::get<ProcessLib::ThermoRichardsMechanics::
355 ConstitutiveStress_StrainTemperature::
356 SwellingDataStateful<DisplacementDim>>(
357 this->current_states_[ip])
360 std::get<
PrevState<ProcessLib::ConstitutiveRelations::
361 MechanicalStrainData<DisplacementDim>>>(
362 this->prev_states_[ip])
365 eps_m_prev.noalias() =
366 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
367 ? eps + C_el.inverse() * sigma_sw
376 DisplacementDim>::assemble(
double const t,
double const dt,
377 std::vector<double>
const& local_x,
378 std::vector<double>
const& local_x_prev,
379 std::vector<double>& local_M_data,
380 std::vector<double>& local_K_data,
381 std::vector<double>& local_rhs_data)
383 assert(local_x.size() == pressure_size + displacement_size);
385 auto const [p_L, u] = localDOF(local_x);
386 auto const [p_L_prev, u_prev] = localDOF(local_x_prev);
389 typename ShapeMatricesTypeDisplacement::template MatrixType<
390 displacement_size + pressure_size,
391 displacement_size + pressure_size>>(
392 local_K_data, displacement_size + pressure_size,
393 displacement_size + pressure_size);
396 typename ShapeMatricesTypeDisplacement::template MatrixType<
397 displacement_size + pressure_size,
398 displacement_size + pressure_size>>(
399 local_M_data, displacement_size + pressure_size,
400 displacement_size + pressure_size);
403 typename ShapeMatricesTypeDisplacement::template VectorType<
404 displacement_size + pressure_size>>(
405 local_rhs_data, displacement_size + pressure_size);
409 DisplacementDim)>::identity2;
412 this->process_data_.media_map.getMedium(this->element_.getID());
413 auto const& liquid_phase = medium->phase(
"AqueousLiquid");
414 auto const& solid_phase = medium->phase(
"Solid");
421 unsigned const n_integration_points =
422 this->integration_method_.getNumberOfPoints();
423 for (
unsigned ip = 0; ip < n_integration_points; ip++)
425 auto const& w = ip_data_[ip].integration_weight;
427 auto const& N_u = ip_data_[ip].N_u;
428 auto const& dNdx_u = ip_data_[ip].dNdx_u;
430 auto const& N_p = ip_data_[ip].N_p;
431 auto const& dNdx_p = ip_data_[ip].dNdx_p;
434 std::nullopt, this->element_.getID(),
438 this->element_, N_u))};
443 ShapeFunctionDisplacement::NPOINTS,
445 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
448 std::get<StrainData<DisplacementDim>>(this->current_states_[ip]);
449 eps.eps.noalias() = B * u;
452 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
453 this->current_states_[ip])
455 auto const S_L_prev =
458 this->prev_states_[ip])
464 double p_cap_prev_ip;
473 auto const temperature =
474 medium->property(MPL::PropertyType::reference_temperature)
475 .template value<double>(variables, x_position, t, dt);
479 medium->property(MPL::PropertyType::biot_coefficient)
480 .template value<double>(variables, x_position, t, dt);
481 auto& SD = this->current_states_[ip];
484 DisplacementDim>>(SD)
490 auto const C_el = ip_data_[ip].computeElasticTangentStiffness(
491 variables, t, x_position, dt, this->solid_material_,
492 *this->material_states_[ip].material_state_variables);
494 auto const beta_SR = (1 - alpha) / this->solid_material_.getBulkModulus(
495 t, x_position, &C_el);
499 liquid_phase.property(MPL::PropertyType::density)
500 .template value<double>(variables, x_position, t, dt);
503 auto const& b = this->process_data_.specific_body_force;
505 S_L = medium->property(MPL::PropertyType::saturation)
506 .template value<double>(variables, x_position, t, dt);
511 double const dS_L_dp_cap =
512 medium->property(MPL::PropertyType::saturation)
513 .template dValue<double>(variables,
514 MPL::Variable::capillary_pressure,
518 double const DeltaS_L_Deltap_cap =
519 (p_cap_ip == p_cap_prev_ip)
521 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
523 auto const chi = [medium, x_position, t, dt](
double const S_L)
527 return medium->property(MPL::PropertyType::bishops_effective_stress)
528 .template value<double>(vs, x_position, t, dt);
530 double const chi_S_L = chi(S_L);
531 double const chi_S_L_prev = chi(S_L_prev);
533 double const p_FR = -chi_S_L * p_cap_ip;
541 auto& phi = std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
542 this->current_states_[ip])
545 auto const phi_prev = std::get<
PrevState<
547 this->prev_states_[ip])
550 phi = medium->property(MPL::PropertyType::porosity)
551 .template value<double>(variables, variables_prev,
559 "RichardsMechanics: Biot-coefficient {} is smaller than "
560 "porosity {} in element/integration point {}/{}.",
561 alpha, phi, this->element_.getID(), ip);
567 std::get<ProcessLib::ThermoRichardsMechanics::
568 ConstitutiveStress_StrainTemperature::
569 SwellingDataStateful<DisplacementDim>>(
570 this->current_states_[ip])
572 auto const& sigma_sw_prev = std::get<
PrevState<
573 ProcessLib::ThermoRichardsMechanics::
574 ConstitutiveStress_StrainTemperature::SwellingDataStateful<
575 DisplacementDim>>>(this->prev_states_[ip])
580 sigma_sw = sigma_sw_prev;
581 if (solid_phase.hasProperty(
582 MPL::PropertyType::swelling_stress_rate))
584 auto const sigma_sw_dot =
587 solid_phase[MPL::PropertyType::swelling_stress_rate]
588 .value(variables, variables_prev, x_position, t,
590 sigma_sw += sigma_sw_dot * dt;
594 identity2.transpose() * C_el.inverse() * sigma_sw;
597 identity2.transpose() * C_el.inverse() * sigma_sw_prev;
607 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
609 auto& transport_porosity =
610 std::get<ProcessLib::ThermoRichardsMechanics::
611 TransportPorosityData>(
612 this->current_states_[ip])
614 auto const transport_porosity_prev =
615 std::get<
PrevState<ProcessLib::ThermoRichardsMechanics::
616 TransportPorosityData>>(
617 this->prev_states_[ip])
622 medium->property(MPL::PropertyType::transport_porosity)
623 .template value<double>(variables, variables_prev,
634 medium->property(MPL::PropertyType::relative_permeability)
635 .template value<double>(variables, x_position, t, dt);
637 liquid_phase.property(MPL::PropertyType::viscosity)
638 .template value<double>(variables, x_position, t, dt);
640 auto const& sigma_sw =
641 std::get<ProcessLib::ThermoRichardsMechanics::
642 ConstitutiveStress_StrainTemperature::
643 SwellingDataStateful<DisplacementDim>>(
644 this->current_states_[ip])
646 auto const& sigma_eff =
648 DisplacementDim>>(this->current_states_[ip])
654 auto const sigma_total =
655 (sigma_eff - alpha * p_FR * identity2).eval();
664 this->material_states_[ip]
665 .material_state_variables->getEquivalentPlasticStrain();
668 medium->property(MPL::PropertyType::permeability)
669 .value(variables, x_position, t, dt));
672 K_intrinsic * rho_LR * k_rel / mu;
678 auto& eps_m = std::get<ProcessLib::ConstitutiveRelations::
679 MechanicalStrainData<DisplacementDim>>(
680 this->current_states_[ip])
683 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
684 ? eps.eps + C_el.inverse() * sigma_sw
692 auto& SD = this->current_states_[ip];
693 auto const& SD_prev = this->prev_states_[ip];
696 DisplacementDim>>(SD);
697 auto const& sigma_eff_prev =
698 std::get<
PrevState<ProcessLib::ConstitutiveRelations::
699 EffectiveStressData<DisplacementDim>>>(
702 std::get<ProcessLib::ConstitutiveRelations::
703 MechanicalStrainData<DisplacementDim>>(SD);
705 std::get<
PrevState<ProcessLib::ConstitutiveRelations::
706 MechanicalStrainData<DisplacementDim>>>(
709 ip_data_[ip].updateConstitutiveRelation(
710 variables, t, x_position, dt, temperature, sigma_eff,
711 sigma_eff_prev, eps_m, eps_m_prev, this->solid_material_,
712 this->material_states_[ip].material_state_variables);
717 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
719 solid_phase.property(MPL::PropertyType::density)
720 .template value<double>(variables, x_position, t, dt);
725 double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
726 rhs.template segment<displacement_size>(displacement_index).noalias() -=
727 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
734 liquid_phase.property(MPL::PropertyType::density)
735 .template dValue<double>(variables,
736 MPL::Variable::liquid_phase_pressure,
739 double const a0 = S_L * (alpha - phi) * beta_SR;
741 double const specific_storage =
742 DeltaS_L_Deltap_cap * (p_cap_ip * a0 - phi) +
743 S_L * (phi * beta_LR + a0);
744 M.template block<pressure_size, pressure_size>(pressure_index,
746 .noalias() += N_p.transpose() * rho_LR * specific_storage * N_p * w;
748 K.template block<pressure_size, pressure_size>(pressure_index,
750 .noalias() += dNdx_p.transpose() * rho_K_over_mu * dNdx_p * w;
752 rhs.template segment<pressure_size>(pressure_index).noalias() +=
753 dNdx_p.transpose() * rho_LR * rho_K_over_mu * b * w;
758 K.template block<displacement_size, pressure_size>(displacement_index,
760 .noalias() -= B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
765 M.template block<pressure_size, displacement_size>(pressure_index,
767 .noalias() += N_p.transpose() * S_L * rho_LR * alpha *
768 identity2.transpose() * B * w;
771 if (this->process_data_.apply_mass_lumping)
773 auto Mpp = M.template block<pressure_size, pressure_size>(
774 pressure_index, pressure_index);
775 Mpp = Mpp.colwise().sum().eval().asDiagonal();
782 ShapeFunctionPressure, DisplacementDim>::
783 assembleWithJacobianEvalConstitutiveSetting(
784 double const t,
double const dt,
787 ShapeFunctionPressure,
788 DisplacementDim>::IpData& ip_data,
795 std::optional<MicroPorosityParameters>
const& micro_porosity_parameters,
801 auto const& liquid_phase = medium->
phase(
"AqueousLiquid");
802 auto const& solid_phase = medium->
phase(
"Solid");
806 DisplacementDim)>::identity2;
808 double const temperature = T_data();
809 double const p_cap_ip = p_cap_data.
p_cap;
810 double const p_cap_prev_ip = p_cap_data.
p_cap_prev;
812 auto const& eps = std::get<StrainData<DisplacementDim>>(SD);
814 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(SD).S_L;
815 auto const S_L_prev =
821 medium->
property(MPL::PropertyType::biot_coefficient)
822 .template value<double>(variables, x_position, t, dt);
823 *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD) = alpha;
827 DisplacementDim>>(SD)
833 auto const C_el = ip_data.computeElasticTangentStiffness(
834 variables, t, x_position, dt, solid_material,
838 (1 - alpha) / solid_material.
getBulkModulus(t, x_position, &C_el);
840 std::get<ProcessLib::ThermoRichardsMechanics::SolidCompressibilityData>(CD)
844 liquid_phase.property(MPL::PropertyType::density)
845 .template value<double>(variables, x_position, t, dt);
847 *std::get<LiquidDensity>(CD) = rho_LR;
849 S_L = medium->
property(MPL::PropertyType::saturation)
850 .template value<double>(variables, x_position, t, dt);
855 double const dS_L_dp_cap =
856 medium->
property(MPL::PropertyType::saturation)
857 .template dValue<double>(variables,
858 MPL::Variable::capillary_pressure,
860 std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(CD)
861 .dS_L_dp_cap = dS_L_dp_cap;
864 double const DeltaS_L_Deltap_cap =
865 (p_cap_ip == p_cap_prev_ip)
867 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
868 std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap =
871 auto const chi = [medium, x_position, t, dt](
double const S_L)
875 return medium->
property(MPL::PropertyType::bishops_effective_stress)
876 .template value<double>(vs, x_position, t, dt);
878 double const chi_S_L = chi(S_L);
879 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD).chi_S_L =
881 double const chi_S_L_prev = chi(S_L_prev);
882 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::BishopsData>>(CD)
883 ->chi_S_L = chi_S_L_prev;
885 auto const dchi_dS_L =
886 medium->
property(MPL::PropertyType::bishops_effective_stress)
887 .template dValue<double>(
888 variables, MPL::Variable::liquid_saturation, x_position, t, dt);
889 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD).dchi_dS_L =
892 double const p_FR = -chi_S_L * p_cap_ip;
904 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(SD).phi;
906 auto const phi_prev =
912 phi = medium->
property(MPL::PropertyType::porosity)
913 .template value<double>(variables, variables_prev, x_position,
917 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi = phi;
923 ?
static_cast<std::ptrdiff_t
>(*x_position.
getElementID())
924 :
static_cast<std::ptrdiff_t
>(-1);
926 "RichardsMechanics: Biot-coefficient {} is smaller than porosity "
931 auto const mu = liquid_phase.property(MPL::PropertyType::viscosity)
932 .template value<double>(variables, x_position, t, dt);
933 *std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(CD) =
939 std::get<ProcessLib::ThermoRichardsMechanics::
940 ConstitutiveStress_StrainTemperature::
941 SwellingDataStateful<DisplacementDim>>(SD);
942 auto const& sigma_sw_prev =
943 std::get<
PrevState<ProcessLib::ThermoRichardsMechanics::
944 ConstitutiveStress_StrainTemperature::
945 SwellingDataStateful<DisplacementDim>>>(
947 auto const transport_porosity_prev = std::get<
PrevState<
950 auto const phi_prev = std::get<
953 auto& transport_porosity = std::get<
955 auto& p_L_m = std::get<MicroPressure>(SD);
956 auto const p_L_m_prev = std::get<PrevState<MicroPressure>>(SD_prev);
957 auto& S_L_m = std::get<MicroSaturation>(SD);
958 auto const S_L_m_prev = std::get<PrevState<MicroSaturation>>(SD_prev);
961 *medium, solid_phase, C_el, rho_LR, mu, micro_porosity_parameters,
962 alpha, phi, p_cap_ip, variables, variables_prev, x_position, t, dt,
963 sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev,
964 transport_porosity, p_L_m_prev, S_L_m_prev, p_L_m, S_L_m);
967 if (medium->
hasProperty(MPL::PropertyType::transport_porosity))
969 if (!medium->
hasProperty(MPL::PropertyType::saturation_micro))
971 auto& transport_porosity =
976 auto const transport_porosity_prev = std::get<
PrevState<
983 medium->
property(MPL::PropertyType::transport_porosity)
984 .template value<double>(variables, variables_prev,
998 auto const sigma_total =
1000 DisplacementDim>>(SD)
1002 alpha * p_FR * identity2)
1011 ->getEquivalentPlasticStrain();
1013 double const k_rel =
1014 medium->
property(MPL::PropertyType::relative_permeability)
1015 .template value<double>(variables, x_position, t, dt);
1018 medium->
property(MPL::PropertyType::permeability)
1019 .
value(variables, x_position, t, dt));
1036 std::get<ProcessLib::ThermoRichardsMechanics::
1037 ConstitutiveStress_StrainTemperature::
1038 SwellingDataStateful<DisplacementDim>>(SD)
1043 DisplacementDim>>(SD)
1046 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
1047 ? eps.eps + C_el.inverse() * sigma_sw
1057 DisplacementDim>>(SD);
1058 auto const& sigma_eff_prev =
1059 std::get<
PrevState<ProcessLib::ConstitutiveRelations::
1060 EffectiveStressData<DisplacementDim>>>(
1064 DisplacementDim>>(SD);
1066 std::get<
PrevState<ProcessLib::ConstitutiveRelations::
1067 MechanicalStrainData<DisplacementDim>>>(
1070 auto C = ip_data.updateConstitutiveRelation(
1071 variables, t, x_position, dt, temperature, sigma_eff,
1072 sigma_eff_prev, eps_m, eps_m_prev, solid_material,
1075 *std::get<StiffnessTensor<DisplacementDim>>(CD) = std::move(C);
1081 DisplacementDim>>(SD)
1082 .sigma_eff.dot(identity2) /
1085 solid_phase.property(MPL::PropertyType::density)
1086 .template value<double>(variables, x_position, t, dt);
1088 double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
1089 *std::get<Density>(CD) = rho;
1095 ShapeFunctionPressure, DisplacementDim>::
1096 assembleWithJacobian(
double const t,
double const dt,
1097 std::vector<double>
const& local_x,
1098 std::vector<double>
const& local_x_prev,
1099 std::vector<double>& local_rhs_data,
1100 std::vector<double>& local_Jac_data)
1102 assert(local_x.size() == pressure_size + displacement_size);
1104 auto const [p_L, u] = localDOF(local_x);
1105 auto const [p_L_prev, u_prev] = localDOF(local_x_prev);
1108 typename ShapeMatricesTypeDisplacement::template MatrixType<
1109 displacement_size + pressure_size,
1110 displacement_size + pressure_size>>(
1111 local_Jac_data, displacement_size + pressure_size,
1112 displacement_size + pressure_size);
1115 typename ShapeMatricesTypeDisplacement::template VectorType<
1116 displacement_size + pressure_size>>(
1117 local_rhs_data, displacement_size + pressure_size);
1121 DisplacementDim)>::identity2;
1124 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1128 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1132 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1136 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1139 typename ShapeMatricesTypeDisplacement::template MatrixType<
1140 displacement_size, pressure_size>
1141 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
1142 displacement_size, pressure_size>::Zero(displacement_size,
1145 typename ShapeMatricesTypeDisplacement::template MatrixType<
1146 pressure_size, displacement_size>
1147 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
1148 pressure_size, displacement_size>::Zero(pressure_size,
1151 auto const& medium =
1152 this->process_data_.media_map.getMedium(this->element_.getID());
1153 auto const& liquid_phase = medium->phase(
"AqueousLiquid");
1154 auto const& solid_phase = medium->phase(
"Solid");
1158 unsigned const n_integration_points =
1159 this->integration_method_.getNumberOfPoints();
1160 for (
unsigned ip = 0; ip < n_integration_points; ip++)
1163 auto& SD = this->current_states_[ip];
1164 auto const& SD_prev = this->prev_states_[ip];
1166 this->process_data_, this->solid_material_);
1168 auto const& w = ip_data_[ip].integration_weight;
1170 auto const& N_u = ip_data_[ip].N_u;
1171 auto const& dNdx_u = ip_data_[ip].dNdx_u;
1173 auto const& N_p = ip_data_[ip].N_p;
1174 auto const& dNdx_p = ip_data_[ip].dNdx_p;
1177 std::nullopt, this->element_.getID(),
1181 this->element_, N_u))};
1186 ShapeFunctionDisplacement::NPOINTS,
1188 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
1193 double p_cap_prev_ip;
1202 auto const temperature =
1203 medium->property(MPL::PropertyType::reference_temperature)
1204 .template value<double>(variables, x_position, t, dt);
1207 std::get<StrainData<DisplacementDim>>(SD).eps.noalias() = B * u;
1209 assembleWithJacobianEvalConstitutiveSetting(
1210 t, dt, x_position, ip_data_[ip], variables, variables_prev, medium,
1213 p_cap_ip, p_cap_prev_ip,
1214 Eigen::Vector<double, DisplacementDim>::Zero()},
1215 CD, SD, SD_prev, this->process_data_.micro_porosity_parameters,
1216 this->solid_material_, this->material_states_[ip]);
1219 auto const& C = *std::get<StiffnessTensor<DisplacementDim>>(CD);
1221 .template block<displacement_size, displacement_size>(
1222 displacement_index, displacement_index)
1223 .noalias() += B.transpose() * C * B * w;
1226 auto const& b = this->process_data_.specific_body_force;
1229 auto const& sigma_eff =
1231 DisplacementDim>>(this->current_states_[ip])
1233 double const rho = *std::get<Density>(CD);
1234 local_rhs.template segment<displacement_size>(displacement_index)
1235 .noalias() -= (B.transpose() * sigma_eff -
1236 N_u_op(N_u).transpose() * rho * b) *
1244 double const alpha =
1245 *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD);
1246 double const dS_L_dp_cap =
1247 std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(
1252 double const chi_S_L =
1253 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1256 B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
1257 double const dchi_dS_L =
1258 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1262 .template block<displacement_size, pressure_size>(
1263 displacement_index, pressure_index)
1264 .noalias() -= B.transpose() * alpha *
1265 (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
1266 identity2 * N_p * w;
1270 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi;
1271 double const rho_LR = *std::get<LiquidDensity>(CD);
1273 .template block<displacement_size, pressure_size>(
1274 displacement_index, pressure_index)
1276 N_u_op(N_u).transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
1287 if (!medium->hasProperty(MPL::PropertyType::saturation_micro) &&
1288 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
1290 using DimMatrix = Eigen::Matrix<double, 3, 3>;
1291 auto const dsigma_sw_dS_L =
1294 .property(MPL::PropertyType::swelling_stress_rate)
1295 .
template dValue<DimMatrix>(
1296 variables, variables_prev,
1297 MPL::Variable::liquid_saturation, x_position, t,
1300 .template block<displacement_size, pressure_size>(
1301 displacement_index, pressure_index)
1303 B.transpose() * dsigma_sw_dS_L * dS_L_dp_cap * N_p * w;
1309 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
1310 this->current_states_[ip])
1312 if (this->process_data_.explicit_hm_coupling_in_unsaturated_zone)
1314 double const chi_S_L_prev = std::get<
PrevState<
1317 Kpu.noalias() += N_p.transpose() * chi_S_L_prev * rho_LR * alpha *
1318 identity2.transpose() * B * w;
1322 Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
1323 identity2.transpose() * B * w;
1330 double const k_rel =
1332 DisplacementDim>>(CD)
1334 auto const& K_intrinsic =
1336 DisplacementDim>>(CD)
1339 *std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(
1344 laplace_p.noalias() +=
1345 dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
1347 auto const beta_LR =
1349 liquid_phase.property(MPL::PropertyType::density)
1350 .template dValue<double>(variables,
1351 MPL::Variable::liquid_phase_pressure,
1354 double const beta_SR =
1359 double const a0 = (alpha - phi) * beta_SR;
1360 double const specific_storage_a_p = S_L * (phi * beta_LR + S_L * a0);
1361 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
1363 double const dspecific_storage_a_p_dp_cap =
1364 dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0);
1365 double const dspecific_storage_a_S_dp_cap =
1366 -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
1368 storage_p_a_p.noalias() +=
1369 N_p.transpose() * rho_LR * specific_storage_a_p * N_p * w;
1371 double const DeltaS_L_Deltap_cap =
1372 std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap;
1373 storage_p_a_S.noalias() -= N_p.transpose() * rho_LR *
1374 specific_storage_a_S * DeltaS_L_Deltap_cap *
1378 .template block<pressure_size, pressure_size>(pressure_index,
1380 .noalias() += N_p.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
1381 rho_LR * dspecific_storage_a_p_dp_cap * N_p * w;
1383 double const S_L_prev =
1386 this->prev_states_[ip])
1388 storage_p_a_S_Jpp.noalias() -=
1389 N_p.transpose() * rho_LR *
1390 ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
1391 specific_storage_a_S * dS_L_dp_cap) /
1394 if (!this->process_data_.explicit_hm_coupling_in_unsaturated_zone)
1397 .template block<pressure_size, pressure_size>(pressure_index,
1399 .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap * alpha *
1400 identity2.transpose() * B * (u - u_prev) / dt *
1404 double const dk_rel_dS_l =
1405 medium->property(MPL::PropertyType::relative_permeability)
1406 .template dValue<double>(variables,
1407 MPL::Variable::liquid_saturation,
1410 grad_p_cap = -dNdx_p * p_L;
1412 .template block<pressure_size, pressure_size>(pressure_index,
1414 .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
1415 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1418 .template block<pressure_size, pressure_size>(pressure_index,
1420 .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
1421 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1423 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
1424 dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
1426 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1428 double const alpha_bar =
1429 this->process_data_.micro_porosity_parameters
1430 ->mass_exchange_coefficient;
1432 *std::get<MicroPressure>(this->current_states_[ip]);
1433 local_rhs.template segment<pressure_size>(pressure_index)
1435 N_p.transpose() * alpha_bar / mu * (-p_cap_ip - p_L_m) * w;
1438 .template block<pressure_size, pressure_size>(pressure_index,
1440 .noalias() += N_p.transpose() * alpha_bar / mu * N_p * w;
1441 if (p_cap_ip != p_cap_prev_ip)
1443 auto const p_L_m_prev = **std::get<PrevState<MicroPressure>>(
1444 this->prev_states_[ip]);
1446 .template block<pressure_size, pressure_size>(
1447 pressure_index, pressure_index)
1448 .noalias() += N_p.transpose() * alpha_bar / mu *
1449 (p_L_m - p_L_m_prev) /
1450 (p_cap_ip - p_cap_prev_ip) * N_p * w;
1455 if (this->process_data_.apply_mass_lumping)
1457 storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
1458 storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
1460 storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
1465 .template block<pressure_size, pressure_size>(pressure_index,
1467 .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
1471 .template block<pressure_size, displacement_size>(pressure_index,
1473 .noalias() = Kpu / dt;
1476 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
1478 (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
1479 Kpu * (u - u_prev) / dt;
1482 local_rhs.template segment<displacement_size>(displacement_index)
1483 .noalias() += Kup * p_L;
1541 ShapeFunctionPressure, DisplacementDim>::
1542 computeSecondaryVariableConcrete(
double const t,
double const dt,
1543 Eigen::VectorXd
const& local_x,
1544 Eigen::VectorXd
const& local_x_prev)
1546 auto const [p_L, u] = localDOF(local_x);
1547 auto const [p_L_prev, u_prev] = localDOF(local_x_prev);
1551 DisplacementDim)>::identity2;
1553 auto const& medium =
1554 this->process_data_.media_map.getMedium(this->element_.getID());
1555 auto const& liquid_phase = medium->phase(
"AqueousLiquid");
1556 auto const& solid_phase = medium->phase(
"Solid");
1560 unsigned const n_integration_points =
1561 this->integration_method_.getNumberOfPoints();
1563 double saturation_avg = 0;
1564 double porosity_avg = 0;
1567 KV sigma_avg = KV::Zero();
1569 for (
unsigned ip = 0; ip < n_integration_points; ip++)
1571 auto const& N_p = ip_data_[ip].N_p;
1572 auto const& N_u = ip_data_[ip].N_u;
1573 auto const& dNdx_u = ip_data_[ip].dNdx_u;
1576 std::nullopt, this->element_.getID(),
1580 this->element_, N_u))};
1585 ShapeFunctionDisplacement::NPOINTS,
1587 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
1592 double p_cap_prev_ip;
1601 auto const temperature =
1602 medium->property(MPL::PropertyType::reference_temperature)
1603 .template value<double>(variables, x_position, t, dt);
1607 std::get<StrainData<DisplacementDim>>(this->current_states_[ip])
1609 eps.noalias() = B * u;
1611 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
1612 this->current_states_[ip])
1614 auto const S_L_prev =
1617 this->prev_states_[ip])
1619 S_L = medium->property(MPL::PropertyType::saturation)
1620 .template value<double>(variables, x_position, t, dt);
1624 auto const chi = [medium, x_position, t, dt](
double const S_L)
1628 return medium->property(MPL::PropertyType::bishops_effective_stress)
1629 .template value<double>(vs, x_position, t, dt);
1631 double const chi_S_L = chi(S_L);
1632 double const chi_S_L_prev = chi(S_L_prev);
1635 medium->property(MPL::PropertyType::biot_coefficient)
1636 .template value<double>(variables, x_position, t, dt);
1637 auto& SD = this->current_states_[ip];
1640 DisplacementDim>>(SD)
1646 auto const C_el = ip_data_[ip].computeElasticTangentStiffness(
1647 variables, t, x_position, dt, this->solid_material_,
1648 *this->material_states_[ip].material_state_variables);
1650 auto const beta_SR = (1 - alpha) / this->solid_material_.getBulkModulus(
1651 t, x_position, &C_el);
1661 auto& phi = std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
1662 this->current_states_[ip])
1665 auto const phi_prev = std::get<
PrevState<
1667 this->prev_states_[ip])
1669 variables_prev.
porosity = phi_prev;
1670 phi = medium->property(MPL::PropertyType::porosity)
1671 .template value<double>(variables, variables_prev,
1677 liquid_phase.property(MPL::PropertyType::density)
1678 .template value<double>(variables, x_position, t, dt);
1681 liquid_phase.property(MPL::PropertyType::viscosity)
1682 .template value<double>(variables, x_position, t, dt);
1687 std::get<ProcessLib::ThermoRichardsMechanics::
1688 ConstitutiveStress_StrainTemperature::
1689 SwellingDataStateful<DisplacementDim>>(
1690 this->current_states_[ip]);
1691 auto const& sigma_sw_prev = std::get<
1692 PrevState<ProcessLib::ThermoRichardsMechanics::
1693 ConstitutiveStress_StrainTemperature::
1694 SwellingDataStateful<DisplacementDim>>>(
1695 this->prev_states_[ip]);
1696 auto const transport_porosity_prev = std::get<
PrevState<
1698 this->prev_states_[ip]);
1699 auto const phi_prev = std::get<
1701 this->prev_states_[ip]);
1702 auto& transport_porosity = std::get<
1704 this->current_states_[ip]);
1705 auto& p_L_m = std::get<MicroPressure>(this->current_states_[ip]);
1706 auto const p_L_m_prev =
1707 std::get<PrevState<MicroPressure>>(this->prev_states_[ip]);
1708 auto& S_L_m = std::get<MicroSaturation>(this->current_states_[ip]);
1709 auto const S_L_m_prev =
1710 std::get<PrevState<MicroSaturation>>(this->prev_states_[ip]);
1713 *medium, solid_phase, C_el, rho_LR, mu,
1714 this->process_data_.micro_porosity_parameters, alpha, phi,
1715 p_cap_ip, variables, variables_prev, x_position, t, dt,
1716 sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev,
1717 transport_porosity, p_L_m_prev, S_L_m_prev, p_L_m, S_L_m);
1720 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
1722 if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
1724 auto& transport_porosity =
1725 std::get<ProcessLib::ThermoRichardsMechanics::
1726 TransportPorosityData>(
1727 this->current_states_[ip])
1729 auto const transport_porosity_prev =
1730 std::get<
PrevState<ProcessLib::ThermoRichardsMechanics::
1731 TransportPorosityData>>(
1732 this->prev_states_[ip])
1737 transport_porosity =
1738 medium->property(MPL::PropertyType::transport_porosity)
1739 .template value<double>(variables, variables_prev,
1749 auto const& sigma_eff =
1751 DisplacementDim>>(this->current_states_[ip])
1757 auto const sigma_total =
1758 (sigma_eff + alpha * chi_S_L * identity2 * p_cap_ip).eval();
1766 this->material_states_[ip]
1767 .material_state_variables->getEquivalentPlasticStrain();
1770 medium->property(MPL::PropertyType::permeability)
1771 .value(variables, x_position, t, dt));
1773 double const k_rel =
1774 medium->property(MPL::PropertyType::relative_permeability)
1775 .template value<double>(variables, x_position, t, dt);
1779 double const p_FR = -chi_S_L * p_cap_ip;
1782 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
1784 solid_phase.property(MPL::PropertyType::density)
1785 .template value<double>(variables, x_position, t, dt);
1786 *std::get<DrySolidDensity>(this->output_data_[ip]) = (1 - phi) * rho_SR;
1789 auto& SD = this->current_states_[ip];
1790 auto const& sigma_sw =
1791 std::get<ProcessLib::ThermoRichardsMechanics::
1792 ConstitutiveStress_StrainTemperature::
1793 SwellingDataStateful<DisplacementDim>>(SD)
1796 std::get<ProcessLib::ConstitutiveRelations::
1797 MechanicalStrainData<DisplacementDim>>(SD)
1800 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
1801 ? eps + C_el.inverse() * sigma_sw
1809 auto& SD = this->current_states_[ip];
1810 auto const& SD_prev = this->prev_states_[ip];
1813 DisplacementDim>>(SD);
1814 auto const& sigma_eff_prev =
1815 std::get<
PrevState<ProcessLib::ConstitutiveRelations::
1816 EffectiveStressData<DisplacementDim>>>(
1819 std::get<ProcessLib::ConstitutiveRelations::
1820 MechanicalStrainData<DisplacementDim>>(SD);
1821 auto const& eps_m_prev =
1822 std::get<
PrevState<ProcessLib::ConstitutiveRelations::
1823 MechanicalStrainData<DisplacementDim>>>(
1826 ip_data_[ip].updateConstitutiveRelation(
1827 variables, t, x_position, dt, temperature, sigma_eff,
1828 sigma_eff_prev, eps_m, eps_m_prev, this->solid_material_,
1829 this->material_states_[ip].material_state_variables);
1832 auto const& b = this->process_data_.specific_body_force;
1835 auto const& dNdx_p = ip_data_[ip].dNdx_p;
1838 this->output_data_[ip])
1839 ->noalias() = -K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b;
1841 saturation_avg += S_L;
1842 porosity_avg += phi;
1843 sigma_avg += sigma_eff;
1845 saturation_avg /= n_integration_points;
1846 porosity_avg /= n_integration_points;
1847 sigma_avg /= n_integration_points;
1849 (*this->process_data_.element_saturation)[this->element_.getID()] =
1851 (*this->process_data_.element_porosity)[this->element_.getID()] =
1855 &(*this->process_data_.element_stresses)[this->element_.getID() *
1856 KV::RowsAtCompileTime]) =
1860 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
1861 DisplacementDim>(this->element_, this->is_axially_symmetric_, p_L,
1862 *this->process_data_.pressure_interpolated);