OGS
RichardsMechanicsFEM-impl.h
Go to the documentation of this file.
1
12#pragma once
13
14#include <Eigen/LU>
15#include <cassert>
16
17#include "ComputeMicroPorosity.h"
18#include "ConstitutiveRelations/ConstitutiveModels.h"
19#include "IntegrationPointData.h"
20#include "MaterialLib/MPL/Medium.h"
21#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
22#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
23#include "MathLib/EigenBlockMatrixView.h"
24#include "MathLib/KelvinVector.h"
25#include "NumLib/Fem/Interpolation.h"
26#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
27#include "ProcessLib/Utils/TransposeInPlace.h"
28#include "RichardsMechanicsFEM.h"
29
30namespace ProcessLib
31{
32namespace RichardsMechanics
33{
34template <int DisplacementDim>
35void updateSwellingStressAndVolumetricStrain(
36 MaterialPropertyLib::Medium const& medium,
37 MaterialPropertyLib::Phase const& solid_phase,
38 MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> const& C_el,
39 double const rho_LR, double const mu,
40 std::optional<MicroPorosityParameters> micro_porosity_parameters,
41 double const alpha, double const phi, double const p_cap_ip,
42 MPL::VariableArray& variables, MPL::VariableArray& variables_prev,
43 ParameterLib::SpatialPosition const& x_position, double const t,
44 double const dt,
45 ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature::
46 SwellingDataStateful<DisplacementDim>& sigma_sw,
47 PrevState<ProcessLib::ThermoRichardsMechanics::
48 ConstitutiveStress_StrainTemperature::SwellingDataStateful<
49 DisplacementDim>> const& sigma_sw_prev,
50 PrevState<ProcessLib::ThermoRichardsMechanics::TransportPorosityData> const
51 phi_M_prev,
52 PrevState<ProcessLib::ThermoRichardsMechanics::PorosityData> const phi_prev,
53 ProcessLib::ThermoRichardsMechanics::TransportPorosityData& phi_M,
54 PrevState<MicroPressure> const p_L_m_prev,
55 PrevState<MicroSaturation> const S_L_m_prev, MicroPressure& p_L_m,
56 MicroSaturation& S_L_m)
57{
58 auto const& identity2 = MathLib::KelvinVector::Invariants<
59 MathLib::KelvinVector::kelvin_vector_dimensions(
60 DisplacementDim)>::identity2;
61
62 if (!medium.hasProperty(MPL::PropertyType::saturation_micro))
63 {
64 // If there is swelling, compute it. Update volumetric strain rate,
65 // s.t. it corresponds to the mechanical part only.
66 sigma_sw = *sigma_sw_prev;
67 if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
68 {
69 auto const sigma_sw_dot =
70 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
71 MPL::formEigenTensor<3>(
72 solid_phase[MPL::PropertyType::swelling_stress_rate]
73 .value(variables, variables_prev, x_position, t,
74 dt)));
75 sigma_sw.sigma_sw += sigma_sw_dot * dt;
76
77 // !!! Misusing volumetric strain for mechanical volumetric
78 // strain just to update the transport porosity !!!
79 variables.volumetric_strain +=
80 identity2.transpose() * C_el.inverse() * sigma_sw.sigma_sw;
81 variables_prev.volumetric_strain += identity2.transpose() *
82 C_el.inverse() *
83 sigma_sw_prev->sigma_sw;
84 }
85 }
86
87 // TODO (naumov) saturation_micro must be always defined together with
88 // the micro_porosity_parameters.
89 if (medium.hasProperty(MPL::PropertyType::saturation_micro))
90 {
91 double const phi_m_prev = phi_prev->phi - phi_M_prev->phi;
92
93 auto const [delta_phi_m, delta_e_sw, delta_p_L_m, delta_sigma_sw] =
94 computeMicroPorosity<DisplacementDim>(
95 identity2.transpose() * C_el.inverse(), rho_LR, mu,
96 *micro_porosity_parameters, alpha, phi, -p_cap_ip, **p_L_m_prev,
97 variables_prev, **S_L_m_prev, phi_m_prev, x_position, t, dt,
98 medium.property(MPL::PropertyType::saturation_micro),
99 solid_phase.property(MPL::PropertyType::swelling_stress_rate));
100
101 phi_M.phi = phi - (phi_m_prev + delta_phi_m);
102 variables_prev.transport_porosity = phi_M_prev->phi;
103 variables.transport_porosity = phi_M.phi;
104
105 *p_L_m = **p_L_m_prev + delta_p_L_m;
106 { // Update micro saturation.
107 MPL::VariableArray variables_prev;
108 variables_prev.capillary_pressure = -**p_L_m_prev;
109 MPL::VariableArray variables;
110 variables.capillary_pressure = -*p_L_m;
111
112 *S_L_m = medium.property(MPL::PropertyType::saturation_micro)
113 .template value<double>(variables, x_position, t, dt);
114 }
115 sigma_sw.sigma_sw = sigma_sw_prev->sigma_sw + delta_sigma_sw;
116 }
117}
118
119template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
120 int DisplacementDim>
121RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
122 ShapeFunctionPressure, DisplacementDim>::
123 RichardsMechanicsLocalAssembler(
124 MeshLib::Element const& e,
125 std::size_t const /*local_matrix_size*/,
126 NumLib::GenericIntegrationMethod const& integration_method,
127 bool const is_axially_symmetric,
128 RichardsMechanicsProcessData<DisplacementDim>& process_data)
129 : LocalAssemblerInterface<DisplacementDim>{
130 e, integration_method, is_axially_symmetric, process_data}
131{
132 unsigned const n_integration_points =
133 this->integration_method_.getNumberOfPoints();
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,
140 ShapeMatricesTypeDisplacement,
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
152 ParameterLib::SpatialPosition x_position;
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 =
160 this->integration_method_.getWeightedPoint(ip).getWeight() *
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<
182 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(
183 this->current_states_[ip])
184 .phi;
185 transport_porosity = porosity;
186 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
187 {
188 transport_porosity =
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}
199
200template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
201 int DisplacementDim>
202void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
203 ShapeFunctionPressure, DisplacementDim>::
204 setInitialConditionsConcrete(Eigen::VectorXd const local_x,
205 double const t,
206 int const /*process_id*/)
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
217 ParameterLib::SpatialPosition x_position;
218 x_position.setElementID(this->element_.getID());
219
220 auto const& solid_phase = medium->phase("Solid");
221
222 unsigned const n_integration_points =
223 this->integration_method_.getNumberOfPoints();
224 for (unsigned ip = 0; ip < n_integration_points; ip++)
225 {
226 x_position.setIntegrationPoint(ip);
227
228 auto const& N_p = _ip_data[ip].N_p;
229
230 double p_cap_ip;
231 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
232
233 variables.capillary_pressure = p_cap_ip;
234 variables.liquid_phase_pressure = -p_cap_ip;
235 // setting pG to 1 atm
236 // TODO : rewrite equations s.t. p_L = pG-p_cap
237 variables.gas_phase_pressure = 1.0e5;
238
239 {
240 auto& p_L_m = std::get<MicroPressure>(this->current_states_[ip]);
241 auto& p_L_m_prev =
242 std::get<PrevState<MicroPressure>>(this->prev_states_[ip]);
243 **p_L_m_prev = -p_cap_ip;
244 *p_L_m = -p_cap_ip;
245 }
246
247 auto const temperature =
248 medium->property(MPL::PropertyType::reference_temperature)
249 .template value<double>(variables, x_position, t, dt);
250 variables.temperature = temperature;
251
252 auto& S_L_prev =
253 std::get<
254 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
255 this->prev_states_[ip])
256 ->S_L;
257 S_L_prev = medium->property(MPL::PropertyType::saturation)
258 .template value<double>(variables, x_position, t, dt);
259
260 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
261 {
262 MPL::VariableArray vars;
263 vars.capillary_pressure = p_cap_ip;
264
265 auto& S_L_m = std::get<MicroSaturation>(this->current_states_[ip]);
266 auto& S_L_m_prev =
267 std::get<PrevState<MicroSaturation>>(this->prev_states_[ip]);
268
269 *S_L_m = medium->property(MPL::PropertyType::saturation_micro)
270 .template value<double>(vars, x_position, t, dt);
271 *S_L_m_prev = S_L_m;
272 }
273
274 // Set eps_m_prev from potentially non-zero eps and sigma_sw from
275 // restart.
276 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
277 t, x_position, dt, temperature, this->solid_material_,
278 *this->material_states_[ip].material_state_variables);
279 auto& eps =
280 std::get<StrainData<DisplacementDim>>(this->current_states_[ip])
281 .eps;
282 auto const& sigma_sw =
283 std::get<ProcessLib::ThermoRichardsMechanics::
284 ConstitutiveStress_StrainTemperature::
285 SwellingDataStateful<DisplacementDim>>(
286 this->current_states_[ip])
287 .sigma_sw;
288 auto& eps_m_prev =
289 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
290 ConstitutiveStress_StrainTemperature::
291 MechanicalStrainData<DisplacementDim>>>(
292 this->prev_states_[ip])
293 ->eps_m;
294
295 eps_m_prev.noalias() =
296 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
297 ? eps + C_el.inverse() * sigma_sw
298 : eps;
299 }
300}
301
302template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
303 int DisplacementDim>
304void RichardsMechanicsLocalAssembler<
305 ShapeFunctionDisplacement, ShapeFunctionPressure,
306 DisplacementDim>::assemble(double const t, double const dt,
307 std::vector<double> const& local_x,
308 std::vector<double> const& local_x_prev,
309 std::vector<double>& local_M_data,
310 std::vector<double>& local_K_data,
311 std::vector<double>& local_rhs_data)
312{
313 assert(local_x.size() == pressure_size + displacement_size);
314
315 auto p_L =
316 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
317 pressure_size> const>(local_x.data() + pressure_index,
318 pressure_size);
319
320 auto u =
321 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
322 displacement_size> const>(local_x.data() + displacement_index,
323 displacement_size);
324
325 auto p_L_prev =
326 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
327 pressure_size> const>(local_x_prev.data() + pressure_index,
328 pressure_size);
329
330 auto u_prev =
331 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
332 displacement_size> const>(local_x_prev.data() + displacement_index,
333 displacement_size);
334
335 auto K = MathLib::createZeroedMatrix<
336 typename ShapeMatricesTypeDisplacement::template MatrixType<
337 displacement_size + pressure_size,
338 displacement_size + pressure_size>>(
339 local_K_data, displacement_size + pressure_size,
340 displacement_size + pressure_size);
341
342 auto M = MathLib::createZeroedMatrix<
343 typename ShapeMatricesTypeDisplacement::template MatrixType<
344 displacement_size + pressure_size,
345 displacement_size + pressure_size>>(
346 local_M_data, displacement_size + pressure_size,
347 displacement_size + pressure_size);
348
349 auto rhs = MathLib::createZeroedVector<
350 typename ShapeMatricesTypeDisplacement::template VectorType<
351 displacement_size + pressure_size>>(
352 local_rhs_data, displacement_size + pressure_size);
353
354 auto const& identity2 = MathLib::KelvinVector::Invariants<
355 MathLib::KelvinVector::kelvin_vector_dimensions(
356 DisplacementDim)>::identity2;
357
358 auto const& medium =
359 this->process_data_.media_map.getMedium(this->element_.getID());
360 auto const& liquid_phase = medium->phase("AqueousLiquid");
361 auto const& solid_phase = medium->phase("Solid");
362 MPL::VariableArray variables;
363 MPL::VariableArray variables_prev;
364
365 ParameterLib::SpatialPosition x_position;
366 x_position.setElementID(this->element_.getID());
367
368 unsigned const n_integration_points =
369 this->integration_method_.getNumberOfPoints();
370 for (unsigned ip = 0; ip < n_integration_points; ip++)
371 {
372 x_position.setIntegrationPoint(ip);
373 auto const& w = _ip_data[ip].integration_weight;
374
375 auto const& N_u = _ip_data[ip].N_u;
376 auto const& dNdx_u = _ip_data[ip].dNdx_u;
377
378 auto const& N_p = _ip_data[ip].N_p;
379 auto const& dNdx_p = _ip_data[ip].dNdx_p;
380
381 auto const x_coord =
382 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
383 ShapeMatricesTypeDisplacement>(
384 this->element_, N_u);
385 auto const B =
386 LinearBMatrix::computeBMatrix<DisplacementDim,
387 ShapeFunctionDisplacement::NPOINTS,
388 typename BMatricesType::BMatrixType>(
389 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
390
391 auto& eps =
392 std::get<StrainData<DisplacementDim>>(this->current_states_[ip]);
393 eps.eps.noalias() = B * u;
394
395 auto& S_L =
396 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
397 this->current_states_[ip])
398 .S_L;
399 auto const S_L_prev =
400 std::get<
401 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
402 this->prev_states_[ip])
403 ->S_L;
404
405 double p_cap_ip;
406 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
407
408 double p_cap_prev_ip;
409 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
410
411 variables.capillary_pressure = p_cap_ip;
412 variables.liquid_phase_pressure = -p_cap_ip;
413 // setting pG to 1 atm
414 // TODO : rewrite equations s.t. p_L = pG-p_cap
415 variables.gas_phase_pressure = 1.0e5;
416
417 auto const temperature =
418 medium->property(MPL::PropertyType::reference_temperature)
419 .template value<double>(variables, x_position, t, dt);
420 variables.temperature = temperature;
421
422 auto const alpha =
423 medium->property(MPL::PropertyType::biot_coefficient)
424 .template value<double>(variables, x_position, t, dt);
425 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
426 t, x_position, dt, temperature, this->solid_material_,
427 *this->material_states_[ip].material_state_variables);
428
429 auto const beta_SR = (1 - alpha) / this->solid_material_.getBulkModulus(
430 t, x_position, &C_el);
431 variables.grain_compressibility = beta_SR;
432
433 auto const rho_LR =
434 liquid_phase.property(MPL::PropertyType::density)
435 .template value<double>(variables, x_position, t, dt);
436 variables.density = rho_LR;
437
438 auto const& b = this->process_data_.specific_body_force;
439
440 S_L = medium->property(MPL::PropertyType::saturation)
441 .template value<double>(variables, x_position, t, dt);
442 variables.liquid_saturation = S_L;
443 variables_prev.liquid_saturation = S_L_prev;
444
445 // tangent derivative for Jacobian
446 double const dS_L_dp_cap =
447 medium->property(MPL::PropertyType::saturation)
448 .template dValue<double>(variables,
449 MPL::Variable::capillary_pressure,
450 x_position, t, dt);
451 // secant derivative from time discretization for storage
452 // use tangent, if secant is not available
453 double const DeltaS_L_Deltap_cap =
454 (p_cap_ip == p_cap_prev_ip)
455 ? dS_L_dp_cap
456 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
457
458 auto const chi = [medium, x_position, t, dt](double const S_L)
459 {
460 MPL::VariableArray vs;
461 vs.liquid_saturation = S_L;
462 return medium->property(MPL::PropertyType::bishops_effective_stress)
463 .template value<double>(vs, x_position, t, dt);
464 };
465 double const chi_S_L = chi(S_L);
466 double const chi_S_L_prev = chi(S_L_prev);
467
468 double const p_FR = -chi_S_L * p_cap_ip;
469 variables.effective_pore_pressure = p_FR;
470 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
471
472 // Set volumetric strain rate for the general case without swelling.
473 variables.volumetric_strain = Invariants::trace(eps.eps);
474 variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
475
476 auto& phi = std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
477 this->current_states_[ip])
478 .phi;
479 { // Porosity update
480 auto const phi_prev = std::get<PrevState<
481 ProcessLib::ThermoRichardsMechanics::PorosityData>>(
482 this->prev_states_[ip])
483 ->phi;
484 variables_prev.porosity = phi_prev;
485 phi = medium->property(MPL::PropertyType::porosity)
486 .template value<double>(variables, variables_prev,
487 x_position, t, dt);
488 variables.porosity = phi;
489 }
490
491 if (alpha < phi)
492 {
493 OGS_FATAL(
494 "RichardsMechanics: Biot-coefficient {} is smaller than "
495 "porosity {} in element/integration point {}/{}.",
496 alpha, phi, this->element_.getID(), ip);
497 }
498
499 // Swelling and possibly volumetric strain rate update.
500 {
501 auto& sigma_sw =
502 std::get<ProcessLib::ThermoRichardsMechanics::
503 ConstitutiveStress_StrainTemperature::
504 SwellingDataStateful<DisplacementDim>>(
505 this->current_states_[ip])
506 .sigma_sw;
507 auto const& sigma_sw_prev = std::get<PrevState<
508 ProcessLib::ThermoRichardsMechanics::
509 ConstitutiveStress_StrainTemperature::SwellingDataStateful<
510 DisplacementDim>>>(this->prev_states_[ip])
511 ->sigma_sw;
512
513 // If there is swelling, compute it. Update volumetric strain rate,
514 // s.t. it corresponds to the mechanical part only.
515 sigma_sw = sigma_sw_prev;
516 if (solid_phase.hasProperty(
517 MPL::PropertyType::swelling_stress_rate))
518 {
519 auto const sigma_sw_dot =
520 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
521 MPL::formEigenTensor<3>(
522 solid_phase[MPL::PropertyType::swelling_stress_rate]
523 .value(variables, variables_prev, x_position, t,
524 dt)));
525 sigma_sw += sigma_sw_dot * dt;
526
527 // !!! Misusing volumetric strain for mechanical volumetric
528 // strain just to update the transport porosity !!!
529 variables.volumetric_strain +=
530 identity2.transpose() * C_el.inverse() * sigma_sw;
531 variables_prev.volumetric_strain +=
532 identity2.transpose() * C_el.inverse() * sigma_sw_prev;
533 }
534
535 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
536 {
537 auto& transport_porosity =
538 std::get<ProcessLib::ThermoRichardsMechanics::
539 TransportPorosityData>(
540 this->current_states_[ip])
541 .phi;
542 auto const transport_porosity_prev =
543 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
544 TransportPorosityData>>(
545 this->prev_states_[ip])
546 ->phi;
547 variables_prev.transport_porosity = transport_porosity_prev;
548
549 transport_porosity =
550 medium->property(MPL::PropertyType::transport_porosity)
551 .template value<double>(variables, variables_prev,
552 x_position, t, dt);
553 variables.transport_porosity = transport_porosity;
554 }
555 else
556 {
557 variables.transport_porosity = phi;
558 }
559 }
560
561 double const k_rel =
562 medium->property(MPL::PropertyType::relative_permeability)
563 .template value<double>(variables, x_position, t, dt);
564 auto const mu =
565 liquid_phase.property(MPL::PropertyType::viscosity)
566 .template value<double>(variables, x_position, t, dt);
567
568 auto const& sigma_sw =
569 std::get<ProcessLib::ThermoRichardsMechanics::
570 ConstitutiveStress_StrainTemperature::
571 SwellingDataStateful<DisplacementDim>>(
572 this->current_states_[ip])
573 .sigma_sw;
574 auto const& sigma_eff =
575 std::get<ProcessLib::ThermoRichardsMechanics::
576 ConstitutiveStress_StrainTemperature::
577 EffectiveStressData<DisplacementDim>>(
578 this->current_states_[ip])
579 .sigma_eff;
580
581 // Set mechanical variables for the intrinsic permeability model
582 // For stress dependent permeability.
583 {
584 auto const sigma_total =
585 (sigma_eff - alpha * p_FR * identity2).eval();
586
587 // For stress dependent permeability.
588 variables.total_stress.emplace<SymmetricTensor>(
589 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
590 sigma_total));
591 }
592
593 variables.equivalent_plastic_strain =
594 this->material_states_[ip]
595 .material_state_variables->getEquivalentPlasticStrain();
596
597 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
598 medium->property(MPL::PropertyType::permeability)
599 .value(variables, x_position, t, dt));
600
601 GlobalDimMatrixType const rho_K_over_mu =
602 K_intrinsic * rho_LR * k_rel / mu;
603
604 //
605 // displacement equation, displacement part
606 //
607 {
608 auto& eps_m =
609 std::get<ProcessLib::ThermoRichardsMechanics::
610 ConstitutiveStress_StrainTemperature::
611 MechanicalStrainData<DisplacementDim>>(
612 this->current_states_[ip])
613 .eps_m;
614 eps_m.noalias() =
615 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
616 ? eps.eps + C_el.inverse() * sigma_sw
617 : eps.eps;
618 variables.mechanical_strain.emplace<
619 MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
620 eps_m);
621 }
622
623 {
624 auto& SD = this->current_states_[ip];
625 auto const& SD_prev = this->prev_states_[ip];
626 auto& sigma_eff =
627 std::get<ProcessLib::ThermoRichardsMechanics::
628 ConstitutiveStress_StrainTemperature::
629 EffectiveStressData<DisplacementDim>>(SD);
630 auto const& sigma_eff_prev = std::get<
631 PrevState<ProcessLib::ThermoRichardsMechanics::
632 ConstitutiveStress_StrainTemperature::
633 EffectiveStressData<DisplacementDim>>>(
634 SD_prev);
635 auto const& eps_m =
636 std::get<ProcessLib::ThermoRichardsMechanics::
637 ConstitutiveStress_StrainTemperature::
638 MechanicalStrainData<DisplacementDim>>(SD);
639 auto& eps_m_prev = std::get<
640 PrevState<ProcessLib::ThermoRichardsMechanics::
641 ConstitutiveStress_StrainTemperature::
642 MechanicalStrainData<DisplacementDim>>>(
643 SD_prev);
644
645 _ip_data[ip].updateConstitutiveRelation(
646 variables, t, x_position, dt, temperature, sigma_eff,
647 sigma_eff_prev, eps_m, eps_m_prev, this->solid_material_,
648 this->material_states_[ip].material_state_variables);
649 }
650
651 // p_SR
652 variables.solid_grain_pressure =
653 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
654 auto const rho_SR =
655 solid_phase.property(MPL::PropertyType::density)
656 .template value<double>(variables, x_position, t, dt);
657
658 //
659 // displacement equation, displacement part
660 //
661 double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
662 rhs.template segment<displacement_size>(displacement_index).noalias() -=
663 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
664
665 //
666 // pressure equation, pressure part.
667 //
668 auto const beta_LR =
669 1 / rho_LR *
670 liquid_phase.property(MPL::PropertyType::density)
671 .template dValue<double>(variables,
672 MPL::Variable::liquid_phase_pressure,
673 x_position, t, dt);
674
675 double const a0 = S_L * (alpha - phi) * beta_SR;
676 // Volumetric average specific storage of the solid and fluid phases.
677 double const specific_storage =
678 DeltaS_L_Deltap_cap * (p_cap_ip * a0 - phi) +
679 S_L * (phi * beta_LR + a0);
680 M.template block<pressure_size, pressure_size>(pressure_index,
681 pressure_index)
682 .noalias() += N_p.transpose() * rho_LR * specific_storage * N_p * w;
683
684 K.template block<pressure_size, pressure_size>(pressure_index,
685 pressure_index)
686 .noalias() += dNdx_p.transpose() * rho_K_over_mu * dNdx_p * w;
687
688 rhs.template segment<pressure_size>(pressure_index).noalias() +=
689 dNdx_p.transpose() * rho_LR * rho_K_over_mu * b * w;
690
691 //
692 // displacement equation, pressure part
693 //
694 K.template block<displacement_size, pressure_size>(displacement_index,
695 pressure_index)
696 .noalias() -= B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
697
698 //
699 // pressure equation, displacement part.
700 //
701 M.template block<pressure_size, displacement_size>(pressure_index,
702 displacement_index)
703 .noalias() += N_p.transpose() * S_L * rho_LR * alpha *
704 identity2.transpose() * B * w;
705 }
706
707 if (this->process_data_.apply_mass_lumping)
708 {
709 auto Mpp = M.template block<pressure_size, pressure_size>(
710 pressure_index, pressure_index);
711 Mpp = Mpp.colwise().sum().eval().asDiagonal();
712 }
713}
714
715template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
716 int DisplacementDim>
717void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
718 ShapeFunctionPressure, DisplacementDim>::
719 assembleWithJacobianEvalConstitutiveSetting(
720 double const t, double const dt,
721 ParameterLib::SpatialPosition const& x_position,
722 RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
723 ShapeFunctionPressure,
724 DisplacementDim>::IpData& ip_data,
725 MPL::VariableArray& variables, MPL::VariableArray& variables_prev,
726 MPL::Medium const* const medium, TemperatureData const T_data,
727 CapillaryPressureData<DisplacementDim> const& p_cap_data,
728 ConstitutiveData<DisplacementDim>& CD,
729 StatefulData<DisplacementDim>& SD,
730 StatefulDataPrev<DisplacementDim> const& SD_prev,
731 std::optional<MicroPorosityParameters> const& micro_porosity_parameters,
732 MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
733 solid_material,
734 ProcessLib::ThermoRichardsMechanics::MaterialStateData<DisplacementDim>&
735 material_state_data)
736{
737 auto const& liquid_phase = medium->phase("AqueousLiquid");
738 auto const& solid_phase = medium->phase("Solid");
739
740 auto const& identity2 = MathLib::KelvinVector::Invariants<
741 MathLib::KelvinVector::kelvin_vector_dimensions(
742 DisplacementDim)>::identity2;
743
744 double const temperature = T_data();
745 double const p_cap_ip = p_cap_data.p_cap;
746 double const p_cap_prev_ip = p_cap_data.p_cap_prev;
747
748 auto const& eps = std::get<StrainData<DisplacementDim>>(SD);
749 auto& S_L =
750 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(SD).S_L;
751 auto const S_L_prev =
752 std::get<
753 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
754 SD_prev)
755 ->S_L;
756 auto const alpha =
757 medium->property(MPL::PropertyType::biot_coefficient)
758 .template value<double>(variables, x_position, t, dt);
759 *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD) = alpha;
760
761 auto const C_el = ip_data.computeElasticTangentStiffness(
762 t, x_position, dt, temperature, solid_material,
763 *material_state_data.material_state_variables);
764
765 auto const beta_SR =
766 (1 - alpha) / solid_material.getBulkModulus(t, x_position, &C_el);
767 variables.grain_compressibility = beta_SR;
768 std::get<ProcessLib::ThermoRichardsMechanics::SolidCompressibilityData>(CD)
769 .beta_SR = beta_SR;
770
771 auto const rho_LR =
772 liquid_phase.property(MPL::PropertyType::density)
773 .template value<double>(variables, x_position, t, dt);
774 variables.density = rho_LR;
775 *std::get<LiquidDensity>(CD) = rho_LR;
776
777 S_L = medium->property(MPL::PropertyType::saturation)
778 .template value<double>(variables, x_position, t, dt);
779 variables.liquid_saturation = S_L;
780 variables_prev.liquid_saturation = S_L_prev;
781
782 // tangent derivative for Jacobian
783 double const dS_L_dp_cap =
784 medium->property(MPL::PropertyType::saturation)
785 .template dValue<double>(variables,
786 MPL::Variable::capillary_pressure,
787 x_position, t, dt);
788 std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(CD)
789 .dS_L_dp_cap = dS_L_dp_cap;
790 // secant derivative from time discretization for storage
791 // use tangent, if secant is not available
792 double const DeltaS_L_Deltap_cap =
793 (p_cap_ip == p_cap_prev_ip)
794 ? dS_L_dp_cap
795 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
796 std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap =
797 DeltaS_L_Deltap_cap;
798
799 auto const chi = [medium, x_position, t, dt](double const S_L)
800 {
801 MPL::VariableArray vs;
802 vs.liquid_saturation = S_L;
803 return medium->property(MPL::PropertyType::bishops_effective_stress)
804 .template value<double>(vs, x_position, t, dt);
805 };
806 double const chi_S_L = chi(S_L);
807 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD).chi_S_L =
808 chi_S_L;
809 double const chi_S_L_prev = chi(S_L_prev);
810 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::BishopsData>>(CD)
811 ->chi_S_L = chi_S_L_prev;
812
813 auto const dchi_dS_L =
814 medium->property(MPL::PropertyType::bishops_effective_stress)
815 .template dValue<double>(
816 variables, MPL::Variable::liquid_saturation, x_position, t, dt);
817 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD).dchi_dS_L =
818 dchi_dS_L;
819
820 double const p_FR = -chi_S_L * p_cap_ip;
821 variables.effective_pore_pressure = p_FR;
822 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
823
824 // Set volumetric strain rate for the general case without swelling.
825 variables.volumetric_strain = Invariants::trace(eps.eps);
826 // TODO (CL) changed that, using eps_prev for the moment, not B * u_prev
827 // variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
828 variables_prev.volumetric_strain = Invariants::trace(
829 std::get<PrevState<StrainData<DisplacementDim>>>(SD_prev)->eps);
830
831 auto& phi =
832 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(SD).phi;
833 { // Porosity update
834 auto const phi_prev =
835 std::get<
836 PrevState<ProcessLib::ThermoRichardsMechanics::PorosityData>>(
837 SD_prev)
838 ->phi;
839 variables_prev.porosity = phi_prev;
840 phi = medium->property(MPL::PropertyType::porosity)
841 .template value<double>(variables, variables_prev, x_position,
842 t, dt);
843 variables.porosity = phi;
844 }
845 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi = phi;
846
847 if (alpha < phi)
848 {
849 auto const eid =
850 x_position.getElementID()
851 ? static_cast<std::ptrdiff_t>(*x_position.getElementID())
852 : static_cast<std::ptrdiff_t>(-1);
853 auto const ip =
854 x_position.getIntegrationPoint()
855 ? static_cast<std::ptrdiff_t>(*x_position.getIntegrationPoint())
856 : static_cast<std::ptrdiff_t>(-1);
857 OGS_FATAL(
858 "RichardsMechanics: Biot-coefficient {} is smaller than "
859 "porosity {} in element/integration point {}/{}.",
860 alpha, phi, eid, ip);
861 }
862
863 auto const mu = liquid_phase.property(MPL::PropertyType::viscosity)
864 .template value<double>(variables, x_position, t, dt);
865 *std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(CD) =
866 mu;
867
868 {
869 // Swelling and possibly volumetric strain rate update.
870 auto& sigma_sw =
871 std::get<ProcessLib::ThermoRichardsMechanics::
872 ConstitutiveStress_StrainTemperature::
873 SwellingDataStateful<DisplacementDim>>(SD);
874 auto const& sigma_sw_prev =
875 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
876 ConstitutiveStress_StrainTemperature::
877 SwellingDataStateful<DisplacementDim>>>(
878 SD_prev);
879 auto const transport_porosity_prev = std::get<PrevState<
880 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>>(
881 SD_prev);
882 auto const phi_prev = std::get<
883 PrevState<ProcessLib::ThermoRichardsMechanics::PorosityData>>(
884 SD_prev);
885 auto& transport_porosity = std::get<
886 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(SD);
887 auto& p_L_m = std::get<MicroPressure>(SD);
888 auto const p_L_m_prev = std::get<PrevState<MicroPressure>>(SD_prev);
889 auto& S_L_m = std::get<MicroSaturation>(SD);
890 auto const S_L_m_prev = std::get<PrevState<MicroSaturation>>(SD_prev);
891
892 updateSwellingStressAndVolumetricStrain<DisplacementDim>(
893 *medium, solid_phase, C_el, rho_LR, mu, micro_porosity_parameters,
894 alpha, phi, p_cap_ip, variables, variables_prev, x_position, t, dt,
895 sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev,
896 transport_porosity, p_L_m_prev, S_L_m_prev, p_L_m, S_L_m);
897 }
898
899 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
900 {
901 if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
902 {
903 auto& transport_porosity =
904 std::get<
905 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(
906 SD)
907 .phi;
908 auto const transport_porosity_prev = std::get<PrevState<
909 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>>(
910 SD_prev)
911 ->phi;
912 variables_prev.transport_porosity = transport_porosity_prev;
913
914 transport_porosity =
915 medium->property(MPL::PropertyType::transport_porosity)
916 .template value<double>(variables, variables_prev,
917 x_position, t, dt);
918 variables.transport_porosity = transport_porosity;
919 }
920 }
921 else
922 {
923 variables.transport_porosity = phi;
924 }
925
926 // Set mechanical variables for the intrinsic permeability model
927 // For stress dependent permeability.
928 {
929 // TODO mechanical constitutive relation will be evaluated afterwards
930 auto const sigma_total =
931 (std::get<ProcessLib::ThermoRichardsMechanics::
932 ConstitutiveStress_StrainTemperature::
933 EffectiveStressData<DisplacementDim>>(SD)
934 .sigma_eff +
935 alpha * p_FR * identity2)
936 .eval();
937 // For stress dependent permeability.
938 variables.total_stress.emplace<SymmetricTensor>(
939 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
940 }
941
942 variables.equivalent_plastic_strain =
943 material_state_data.material_state_variables
944 ->getEquivalentPlasticStrain();
945
946 double const k_rel =
947 medium->property(MPL::PropertyType::relative_permeability)
948 .template value<double>(variables, x_position, t, dt);
949
950 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
951 medium->property(MPL::PropertyType::permeability)
952 .value(variables, x_position, t, dt));
953
954 std::get<
955 ProcessLib::ThermoRichardsMechanics::PermeabilityData<DisplacementDim>>(
956 CD)
957 .k_rel = k_rel;
958 std::get<
959 ProcessLib::ThermoRichardsMechanics::PermeabilityData<DisplacementDim>>(
960 CD)
961 .Ki = K_intrinsic;
962
963 //
964 // displacement equation, displacement part
965 //
966
967 {
968 auto& sigma_sw =
969 std::get<ProcessLib::ThermoRichardsMechanics::
970 ConstitutiveStress_StrainTemperature::
971 SwellingDataStateful<DisplacementDim>>(SD)
972 .sigma_sw;
973
974 auto& eps_m =
975 std::get<ProcessLib::ThermoRichardsMechanics::
976 ConstitutiveStress_StrainTemperature::
977 MechanicalStrainData<DisplacementDim>>(SD)
978 .eps_m;
979 eps_m.noalias() =
980 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
981 ? eps.eps + C_el.inverse() * sigma_sw
982 : eps.eps;
983 variables.mechanical_strain
984 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
985 eps_m);
986 }
987
988 {
989 auto& sigma_eff =
990 std::get<ProcessLib::ThermoRichardsMechanics::
991 ConstitutiveStress_StrainTemperature::
992 EffectiveStressData<DisplacementDim>>(SD);
993 auto const& sigma_eff_prev =
994 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
995 ConstitutiveStress_StrainTemperature::
996 EffectiveStressData<DisplacementDim>>>(
997 SD_prev);
998 auto const& eps_m =
999 std::get<ProcessLib::ThermoRichardsMechanics::
1000 ConstitutiveStress_StrainTemperature::
1001 MechanicalStrainData<DisplacementDim>>(SD);
1002 auto& eps_m_prev =
1003 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
1004 ConstitutiveStress_StrainTemperature::
1005 MechanicalStrainData<DisplacementDim>>>(
1006 SD_prev);
1007
1008 auto C = ip_data.updateConstitutiveRelation(
1009 variables, t, x_position, dt, temperature, sigma_eff,
1010 sigma_eff_prev, eps_m, eps_m_prev, solid_material,
1011 material_state_data.material_state_variables);
1012
1013 *std::get<StiffnessTensor<DisplacementDim>>(CD) = std::move(C);
1014 }
1015
1016 // p_SR
1017 variables.solid_grain_pressure =
1018 p_FR -
1019 std::get<ProcessLib::ThermoRichardsMechanics::
1020 ConstitutiveStress_StrainTemperature::EffectiveStressData<
1021 DisplacementDim>>(SD)
1022 .sigma_eff.dot(identity2) /
1023 (3 * (1 - phi));
1024 auto const rho_SR =
1025 solid_phase.property(MPL::PropertyType::density)
1026 .template value<double>(variables, x_position, t, dt);
1027
1028 double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
1029 *std::get<Density>(CD) = rho;
1030}
1031
1032template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1033 int DisplacementDim>
1034void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
1035 ShapeFunctionPressure, DisplacementDim>::
1036 assembleWithJacobian(double const t, double const dt,
1037 std::vector<double> const& local_x,
1038 std::vector<double> const& local_x_prev,
1039 std::vector<double>& local_rhs_data,
1040 std::vector<double>& local_Jac_data)
1041{
1042 assert(local_x.size() == pressure_size + displacement_size);
1043
1044 auto p_L =
1045 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
1046 pressure_size> const>(local_x.data() + pressure_index,
1047 pressure_size);
1048
1049 auto u =
1050 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
1051 displacement_size> const>(local_x.data() + displacement_index,
1052 displacement_size);
1053
1054 auto p_L_prev =
1055 Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
1056 pressure_size> const>(local_x_prev.data() + pressure_index,
1057 pressure_size);
1058 auto u_prev =
1059 Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
1060 displacement_size> const>(local_x_prev.data() + displacement_index,
1061 displacement_size);
1062
1063 auto local_Jac = MathLib::createZeroedMatrix<
1064 typename ShapeMatricesTypeDisplacement::template MatrixType<
1065 displacement_size + pressure_size,
1066 displacement_size + pressure_size>>(
1067 local_Jac_data, displacement_size + pressure_size,
1068 displacement_size + pressure_size);
1069
1070 auto local_rhs = MathLib::createZeroedVector<
1071 typename ShapeMatricesTypeDisplacement::template VectorType<
1072 displacement_size + pressure_size>>(
1073 local_rhs_data, displacement_size + pressure_size);
1074
1075 auto const& identity2 = MathLib::KelvinVector::Invariants<
1076 MathLib::KelvinVector::kelvin_vector_dimensions(
1077 DisplacementDim)>::identity2;
1078
1079 typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
1080 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1081 pressure_size);
1082
1083 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_p =
1084 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1085 pressure_size);
1086
1087 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S_Jpp =
1088 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1089 pressure_size);
1090
1091 typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S =
1092 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
1093 pressure_size);
1094
1095 typename ShapeMatricesTypeDisplacement::template MatrixType<
1096 displacement_size, pressure_size>
1097 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
1098 displacement_size, pressure_size>::Zero(displacement_size,
1099 pressure_size);
1100
1101 typename ShapeMatricesTypeDisplacement::template MatrixType<
1102 pressure_size, displacement_size>
1103 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
1104 pressure_size, displacement_size>::Zero(pressure_size,
1105 displacement_size);
1106
1107 auto const& medium =
1108 this->process_data_.media_map.getMedium(this->element_.getID());
1109 auto const& liquid_phase = medium->phase("AqueousLiquid");
1110 auto const& solid_phase = medium->phase("Solid");
1111 MPL::VariableArray variables;
1112 MPL::VariableArray variables_prev;
1113
1114 ParameterLib::SpatialPosition x_position;
1115 x_position.setElementID(this->element_.getID());
1116
1117 unsigned const n_integration_points =
1118 this->integration_method_.getNumberOfPoints();
1119 for (unsigned ip = 0; ip < n_integration_points; ip++)
1120 {
1121 ConstitutiveData<DisplacementDim> CD;
1122 auto& SD = this->current_states_[ip];
1123 auto const& SD_prev = this->prev_states_[ip];
1124 [[maybe_unused]] auto models = createConstitutiveModels(
1125 this->process_data_, this->solid_material_);
1126
1127 x_position.setIntegrationPoint(ip);
1128 auto const& w = _ip_data[ip].integration_weight;
1129
1130 auto const& N_u = _ip_data[ip].N_u;
1131 auto const& dNdx_u = _ip_data[ip].dNdx_u;
1132
1133 auto const& N_p = _ip_data[ip].N_p;
1134 auto const& dNdx_p = _ip_data[ip].dNdx_p;
1135
1136 auto const x_coord =
1137 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1138 ShapeMatricesTypeDisplacement>(
1139 this->element_, N_u);
1140 auto const B =
1141 LinearBMatrix::computeBMatrix<DisplacementDim,
1142 ShapeFunctionDisplacement::NPOINTS,
1143 typename BMatricesType::BMatrixType>(
1144 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
1145
1146 double p_cap_ip;
1147 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
1148
1149 double p_cap_prev_ip;
1150 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
1151
1152 variables.capillary_pressure = p_cap_ip;
1153 variables.liquid_phase_pressure = -p_cap_ip;
1154 // setting pG to 1 atm
1155 // TODO : rewrite equations s.t. p_L = pG-p_cap
1156 variables.gas_phase_pressure = 1.0e5;
1157
1158 auto const temperature =
1159 medium->property(MPL::PropertyType::reference_temperature)
1160 .template value<double>(variables, x_position, t, dt);
1161 variables.temperature = temperature;
1162
1163 std::get<StrainData<DisplacementDim>>(SD).eps.noalias() = B * u;
1164
1165 assembleWithJacobianEvalConstitutiveSetting(
1166 t, dt, x_position, _ip_data[ip], variables, variables_prev, medium,
1167 TemperatureData{temperature},
1168 CapillaryPressureData<DisplacementDim>{
1169 p_cap_ip, p_cap_prev_ip,
1170 Eigen::Vector<double, DisplacementDim>::Zero()},
1171 CD, SD, SD_prev, this->process_data_.micro_porosity_parameters,
1172 this->solid_material_, this->material_states_[ip]);
1173
1174 {
1175 auto const& C = *std::get<StiffnessTensor<DisplacementDim>>(CD);
1176 local_Jac
1177 .template block<displacement_size, displacement_size>(
1178 displacement_index, displacement_index)
1179 .noalias() += B.transpose() * C * B * w;
1180 }
1181
1182 auto const& b = this->process_data_.specific_body_force;
1183
1184 {
1185 auto const& sigma_eff =
1186 std::get<ProcessLib::ThermoRichardsMechanics::
1187 ConstitutiveStress_StrainTemperature::
1188 EffectiveStressData<DisplacementDim>>(
1189 this->current_states_[ip])
1190 .sigma_eff;
1191 double const rho = *std::get<Density>(CD);
1192 local_rhs.template segment<displacement_size>(displacement_index)
1193 .noalias() -= (B.transpose() * sigma_eff -
1194 N_u_op(N_u).transpose() * rho * b) *
1195 w;
1196 }
1197
1198 //
1199 // displacement equation, pressure part
1200 //
1201
1202 double const alpha =
1203 *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD);
1204 double const dS_L_dp_cap =
1205 std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(
1206 CD)
1207 .dS_L_dp_cap;
1208
1209 {
1210 double const chi_S_L =
1211 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1212 .chi_S_L;
1213 Kup.noalias() +=
1214 B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
1215 double const dchi_dS_L =
1216 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1217 .dchi_dS_L;
1218
1219 local_Jac
1220 .template block<displacement_size, pressure_size>(
1221 displacement_index, pressure_index)
1222 .noalias() -= B.transpose() * alpha *
1223 (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
1224 identity2 * N_p * w;
1225 }
1226
1227 double const phi =
1228 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi;
1229 double const rho_LR = *std::get<LiquidDensity>(CD);
1230 local_Jac
1231 .template block<displacement_size, pressure_size>(
1232 displacement_index, pressure_index)
1233 .noalias() +=
1234 N_u_op(N_u).transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
1235
1236 // For the swelling stress with double structure model the corresponding
1237 // Jacobian u-p entry would be required, but it does not improve
1238 // convergence and sometimes worsens it:
1239 // if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1240 // {
1241 // -B.transpose() *
1242 // dsigma_sw_dS_L_m* dS_L_m_dp_cap_m*(p_L_m - p_L_m_prev) /
1243 // (p_cap_ip - p_cap_prev_ip) * N_p* w;
1244 // }
1245 if (!medium->hasProperty(MPL::PropertyType::saturation_micro) &&
1246 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
1247 {
1248 using DimMatrix = Eigen::Matrix<double, 3, 3>;
1249 auto const dsigma_sw_dS_L =
1250 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
1251 solid_phase
1252 .property(MPL::PropertyType::swelling_stress_rate)
1253 .template dValue<DimMatrix>(
1254 variables, variables_prev,
1255 MPL::Variable::liquid_saturation, x_position, t,
1256 dt));
1257 local_Jac
1258 .template block<displacement_size, pressure_size>(
1259 displacement_index, pressure_index)
1260 .noalias() +=
1261 B.transpose() * dsigma_sw_dS_L * dS_L_dp_cap * N_p * w;
1262 }
1263 //
1264 // pressure equation, displacement part.
1265 //
1266 double const S_L =
1267 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
1268 this->current_states_[ip])
1269 .S_L;
1270 if (this->process_data_.explicit_hm_coupling_in_unsaturated_zone)
1271 {
1272 double const chi_S_L_prev = std::get<PrevState<
1273 ProcessLib::ThermoRichardsMechanics::BishopsData>>(CD)
1274 ->chi_S_L;
1275 Kpu.noalias() += N_p.transpose() * chi_S_L_prev * rho_LR * alpha *
1276 identity2.transpose() * B * w;
1277 }
1278 else
1279 {
1280 Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
1281 identity2.transpose() * B * w;
1282 }
1283
1284 //
1285 // pressure equation, pressure part.
1286 //
1287
1288 double const k_rel =
1289 std::get<ProcessLib::ThermoRichardsMechanics::PermeabilityData<
1290 DisplacementDim>>(CD)
1291 .k_rel;
1292 auto const& K_intrinsic =
1293 std::get<ProcessLib::ThermoRichardsMechanics::PermeabilityData<
1294 DisplacementDim>>(CD)
1295 .Ki;
1296 double const mu =
1297 *std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(
1298 CD);
1299
1300 GlobalDimMatrixType const rho_Ki_over_mu = K_intrinsic * rho_LR / mu;
1301
1302 laplace_p.noalias() +=
1303 dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
1304
1305 auto const beta_LR =
1306 1 / rho_LR *
1307 liquid_phase.property(MPL::PropertyType::density)
1308 .template dValue<double>(variables,
1309 MPL::Variable::liquid_phase_pressure,
1310 x_position, t, dt);
1311
1312 double const beta_SR =
1313 std::get<
1314 ProcessLib::ThermoRichardsMechanics::SolidCompressibilityData>(
1315 CD)
1316 .beta_SR;
1317 double const a0 = (alpha - phi) * beta_SR;
1318 double const specific_storage_a_p = S_L * (phi * beta_LR + S_L * a0);
1319 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
1320
1321 double const dspecific_storage_a_p_dp_cap =
1322 dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0);
1323 double const dspecific_storage_a_S_dp_cap =
1324 -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
1325
1326 storage_p_a_p.noalias() +=
1327 N_p.transpose() * rho_LR * specific_storage_a_p * N_p * w;
1328
1329 double const DeltaS_L_Deltap_cap =
1330 std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap;
1331 storage_p_a_S.noalias() -= N_p.transpose() * rho_LR *
1332 specific_storage_a_S * DeltaS_L_Deltap_cap *
1333 N_p * w;
1334
1335 local_Jac
1336 .template block<pressure_size, pressure_size>(pressure_index,
1337 pressure_index)
1338 .noalias() += N_p.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
1339 rho_LR * dspecific_storage_a_p_dp_cap * N_p * w;
1340
1341 double const S_L_prev =
1342 std::get<
1343 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
1344 this->prev_states_[ip])
1345 ->S_L;
1346 storage_p_a_S_Jpp.noalias() -=
1347 N_p.transpose() * rho_LR *
1348 ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
1349 specific_storage_a_S * dS_L_dp_cap) /
1350 dt * N_p * w;
1351
1352 if (!this->process_data_.explicit_hm_coupling_in_unsaturated_zone)
1353 {
1354 local_Jac
1355 .template block<pressure_size, pressure_size>(pressure_index,
1356 pressure_index)
1357 .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap * alpha *
1358 identity2.transpose() * B * (u - u_prev) / dt *
1359 N_p * w;
1360 }
1361
1362 double const dk_rel_dS_l =
1363 medium->property(MPL::PropertyType::relative_permeability)
1364 .template dValue<double>(variables,
1365 MPL::Variable::liquid_saturation,
1366 x_position, t, dt);
1367 typename ShapeMatricesTypeDisplacement::GlobalDimVectorType const
1368 grad_p_cap = -dNdx_p * p_L;
1369 local_Jac
1370 .template block<pressure_size, pressure_size>(pressure_index,
1371 pressure_index)
1372 .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
1373 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1374
1375 local_Jac
1376 .template block<pressure_size, pressure_size>(pressure_index,
1377 pressure_index)
1378 .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
1379 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1380
1381 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
1382 dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
1383
1384 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1385 {
1386 double const alpha_bar =
1387 this->process_data_.micro_porosity_parameters
1388 ->mass_exchange_coefficient;
1389 auto const p_L_m =
1390 *std::get<MicroPressure>(this->current_states_[ip]);
1391 local_rhs.template segment<pressure_size>(pressure_index)
1392 .noalias() -=
1393 N_p.transpose() * alpha_bar / mu * (-p_cap_ip - p_L_m) * w;
1394
1395 local_Jac
1396 .template block<pressure_size, pressure_size>(pressure_index,
1397 pressure_index)
1398 .noalias() += N_p.transpose() * alpha_bar / mu * N_p * w;
1399 if (p_cap_ip != p_cap_prev_ip)
1400 {
1401 auto const p_L_m_prev = **std::get<PrevState<MicroPressure>>(
1402 this->prev_states_[ip]);
1403 local_Jac
1404 .template block<pressure_size, pressure_size>(
1405 pressure_index, pressure_index)
1406 .noalias() += N_p.transpose() * alpha_bar / mu *
1407 (p_L_m - p_L_m_prev) /
1408 (p_cap_ip - p_cap_prev_ip) * N_p * w;
1409 }
1410 }
1411 }
1412
1413 if (this->process_data_.apply_mass_lumping)
1414 {
1415 storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
1416 storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
1417 storage_p_a_S_Jpp =
1418 storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
1419 }
1420
1421 // pressure equation, pressure part.
1422 local_Jac
1423 .template block<pressure_size, pressure_size>(pressure_index,
1424 pressure_index)
1425 .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
1426
1427 // pressure equation, displacement part.
1428 local_Jac
1429 .template block<pressure_size, displacement_size>(pressure_index,
1430 displacement_index)
1431 .noalias() = Kpu / dt;
1432
1433 // pressure equation
1434 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
1435 laplace_p * p_L +
1436 (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
1437 Kpu * (u - u_prev) / dt;
1438
1439 // displacement equation
1440 local_rhs.template segment<displacement_size>(displacement_index)
1441 .noalias() += Kup * p_L;
1442}
1443
1444template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1445 int DisplacementDim>
1446void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
1447 ShapeFunctionPressure, DisplacementDim>::
1448 assembleWithJacobianForPressureEquations(
1449 const double /*t*/, double const /*dt*/,
1450 Eigen::VectorXd const& /*local_x*/,
1451 Eigen::VectorXd const& /*local_x_prev*/,
1452 std::vector<double>& /*local_b_data*/,
1453 std::vector<double>& /*local_Jac_data*/)
1454{
1455 OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
1456}
1457
1458template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1459 int DisplacementDim>
1460void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
1461 ShapeFunctionPressure, DisplacementDim>::
1462 assembleWithJacobianForDeformationEquations(
1463 const double /*t*/, double const /*dt*/,
1464 Eigen::VectorXd const& /*local_x*/,
1465 Eigen::VectorXd const& /*local_x_prev*/,
1466 std::vector<double>& /*local_b_data*/,
1467 std::vector<double>& /*local_Jac_data*/)
1468{
1469 OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
1470}
1471
1472template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1473 int DisplacementDim>
1474void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
1475 ShapeFunctionPressure, DisplacementDim>::
1476 assembleWithJacobianForStaggeredScheme(double const t, double const dt,
1477 Eigen::VectorXd const& local_x,
1478 Eigen::VectorXd const& local_x_prev,
1479 int const process_id,
1480 std::vector<double>& local_b_data,
1481 std::vector<double>& local_Jac_data)
1482{
1483 // For the equations with pressure
1484 if (process_id == 0)
1485 {
1486 assembleWithJacobianForPressureEquations(t, dt, local_x, local_x_prev,
1487 local_b_data, local_Jac_data);
1488 return;
1489 }
1490
1491 // For the equations with deformation
1492 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_x_prev,
1493 local_b_data, local_Jac_data);
1494}
1495
1496template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1497 int DisplacementDim>
1498void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
1499 ShapeFunctionPressure, DisplacementDim>::
1500 computeSecondaryVariableConcrete(double const t, double const dt,
1501 Eigen::VectorXd const& local_x,
1502 Eigen::VectorXd const& local_x_prev)
1503{
1504 auto p_L = local_x.template segment<pressure_size>(pressure_index);
1505 auto u = local_x.template segment<displacement_size>(displacement_index);
1506
1507 auto p_L_prev =
1508 local_x_prev.template segment<pressure_size>(pressure_index);
1509 auto u_prev =
1510 local_x_prev.template segment<displacement_size>(displacement_index);
1511
1512 auto const& identity2 = MathLib::KelvinVector::Invariants<
1513 MathLib::KelvinVector::kelvin_vector_dimensions(
1514 DisplacementDim)>::identity2;
1515
1516 auto const& medium =
1517 this->process_data_.media_map.getMedium(this->element_.getID());
1518 auto const& liquid_phase = medium->phase("AqueousLiquid");
1519 auto const& solid_phase = medium->phase("Solid");
1520 MPL::VariableArray variables;
1521 MPL::VariableArray variables_prev;
1522
1523 ParameterLib::SpatialPosition x_position;
1524 x_position.setElementID(this->element_.getID());
1525
1526 unsigned const n_integration_points =
1527 this->integration_method_.getNumberOfPoints();
1528
1529 double saturation_avg = 0;
1530 double porosity_avg = 0;
1531
1532 using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
1533 KV sigma_avg = KV::Zero();
1534
1535 for (unsigned ip = 0; ip < n_integration_points; ip++)
1536 {
1537 x_position.setIntegrationPoint(ip);
1538 auto const& N_p = _ip_data[ip].N_p;
1539 auto const& N_u = _ip_data[ip].N_u;
1540 auto const& dNdx_u = _ip_data[ip].dNdx_u;
1541
1542 auto const x_coord =
1543 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1544 ShapeMatricesTypeDisplacement>(
1545 this->element_, N_u);
1546 auto const B =
1547 LinearBMatrix::computeBMatrix<DisplacementDim,
1548 ShapeFunctionDisplacement::NPOINTS,
1549 typename BMatricesType::BMatrixType>(
1550 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
1551
1552 double p_cap_ip;
1553 NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
1554
1555 double p_cap_prev_ip;
1556 NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
1557
1558 variables.capillary_pressure = p_cap_ip;
1559 variables.liquid_phase_pressure = -p_cap_ip;
1560 // setting pG to 1 atm
1561 // TODO : rewrite equations s.t. p_L = pG-p_cap
1562 variables.gas_phase_pressure = 1.0e5;
1563
1564 auto const temperature =
1565 medium->property(MPL::PropertyType::reference_temperature)
1566 .template value<double>(variables, x_position, t, dt);
1567 variables.temperature = temperature;
1568
1569 auto& eps =
1570 std::get<StrainData<DisplacementDim>>(this->current_states_[ip])
1571 .eps;
1572 eps.noalias() = B * u;
1573 auto& S_L =
1574 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
1575 this->current_states_[ip])
1576 .S_L;
1577 auto const S_L_prev =
1578 std::get<
1579 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
1580 this->prev_states_[ip])
1581 ->S_L;
1582 S_L = medium->property(MPL::PropertyType::saturation)
1583 .template value<double>(variables, x_position, t, dt);
1584 variables.liquid_saturation = S_L;
1585 variables_prev.liquid_saturation = S_L_prev;
1586
1587 auto const chi = [medium, x_position, t, dt](double const S_L)
1588 {
1589 MPL::VariableArray vs;
1590 vs.liquid_saturation = S_L;
1591 return medium->property(MPL::PropertyType::bishops_effective_stress)
1592 .template value<double>(vs, x_position, t, dt);
1593 };
1594 double const chi_S_L = chi(S_L);
1595 double const chi_S_L_prev = chi(S_L_prev);
1596
1597 auto const alpha =
1598 medium->property(MPL::PropertyType::biot_coefficient)
1599 .template value<double>(variables, x_position, t, dt);
1600
1601 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
1602 t, x_position, dt, temperature, this->solid_material_,
1603 *this->material_states_[ip].material_state_variables);
1604
1605 auto const beta_SR = (1 - alpha) / this->solid_material_.getBulkModulus(
1606 t, x_position, &C_el);
1607 variables.grain_compressibility = beta_SR;
1608
1609 variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
1610 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
1611
1612 // Set volumetric strain rate for the general case without swelling.
1613 variables.volumetric_strain = Invariants::trace(eps);
1614 variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
1615
1616 auto& phi = std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
1617 this->current_states_[ip])
1618 .phi;
1619 { // Porosity update
1620 auto const phi_prev = std::get<PrevState<
1621 ProcessLib::ThermoRichardsMechanics::PorosityData>>(
1622 this->prev_states_[ip])
1623 ->phi;
1624 variables_prev.porosity = phi_prev;
1625 phi = medium->property(MPL::PropertyType::porosity)
1626 .template value<double>(variables, variables_prev,
1627 x_position, t, dt);
1628 variables.porosity = phi;
1629 }
1630
1631 auto const rho_LR =
1632 liquid_phase.property(MPL::PropertyType::density)
1633 .template value<double>(variables, x_position, t, dt);
1634 variables.density = rho_LR;
1635 auto const mu =
1636 liquid_phase.property(MPL::PropertyType::viscosity)
1637 .template value<double>(variables, x_position, t, dt);
1638
1639 {
1640 // Swelling and possibly volumetric strain rate update.
1641 auto& sigma_sw =
1642 std::get<ProcessLib::ThermoRichardsMechanics::
1643 ConstitutiveStress_StrainTemperature::
1644 SwellingDataStateful<DisplacementDim>>(
1645 this->current_states_[ip]);
1646 auto const& sigma_sw_prev = std::get<
1647 PrevState<ProcessLib::ThermoRichardsMechanics::
1648 ConstitutiveStress_StrainTemperature::
1649 SwellingDataStateful<DisplacementDim>>>(
1650 this->prev_states_[ip]);
1651 auto const transport_porosity_prev = std::get<PrevState<
1652 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>>(
1653 this->prev_states_[ip]);
1654 auto const phi_prev = std::get<
1655 PrevState<ProcessLib::ThermoRichardsMechanics::PorosityData>>(
1656 this->prev_states_[ip]);
1657 auto& transport_porosity = std::get<
1658 ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(
1659 this->current_states_[ip]);
1660 auto& p_L_m = std::get<MicroPressure>(this->current_states_[ip]);
1661 auto const p_L_m_prev =
1662 std::get<PrevState<MicroPressure>>(this->prev_states_[ip]);
1663 auto& S_L_m = std::get<MicroSaturation>(this->current_states_[ip]);
1664 auto const S_L_m_prev =
1665 std::get<PrevState<MicroSaturation>>(this->prev_states_[ip]);
1666
1667 updateSwellingStressAndVolumetricStrain<DisplacementDim>(
1668 *medium, solid_phase, C_el, rho_LR, mu,
1669 this->process_data_.micro_porosity_parameters, alpha, phi,
1670 p_cap_ip, variables, variables_prev, x_position, t, dt,
1671 sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev,
1672 transport_porosity, p_L_m_prev, S_L_m_prev, p_L_m, S_L_m);
1673 }
1674
1675 if (medium->hasProperty(MPL::PropertyType::transport_porosity))
1676 {
1677 if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
1678 {
1679 auto& transport_porosity =
1680 std::get<ProcessLib::ThermoRichardsMechanics::
1681 TransportPorosityData>(
1682 this->current_states_[ip])
1683 .phi;
1684 auto const transport_porosity_prev =
1685 std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
1686 TransportPorosityData>>(
1687 this->prev_states_[ip])
1688 ->phi;
1689
1690 variables_prev.transport_porosity = transport_porosity_prev;
1691
1692 transport_porosity =
1693 medium->property(MPL::PropertyType::transport_porosity)
1694 .template value<double>(variables, variables_prev,
1695 x_position, t, dt);
1696 variables.transport_porosity = transport_porosity;
1697 }
1698 }
1699 else
1700 {
1701 variables.transport_porosity = phi;
1702 }
1703
1704 auto const& sigma_eff =
1705 std::get<ProcessLib::ThermoRichardsMechanics::
1706 ConstitutiveStress_StrainTemperature::
1707 EffectiveStressData<DisplacementDim>>(
1708 this->current_states_[ip])
1709 .sigma_eff;
1710
1711 // Set mechanical variables for the intrinsic permeability model
1712 // For stress dependent permeability.
1713 {
1714 auto const sigma_total =
1715 (sigma_eff + alpha * chi_S_L * identity2 * p_cap_ip).eval();
1716 // For stress dependent permeability.
1717 variables.total_stress.emplace<SymmetricTensor>(
1718 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
1719 sigma_total));
1720 }
1721
1722 variables.equivalent_plastic_strain =
1723 this->material_states_[ip]
1724 .material_state_variables->getEquivalentPlasticStrain();
1725
1726 auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
1727 medium->property(MPL::PropertyType::permeability)
1728 .value(variables, x_position, t, dt));
1729
1730 double const k_rel =
1731 medium->property(MPL::PropertyType::relative_permeability)
1732 .template value<double>(variables, x_position, t, dt);
1733
1734 GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
1735
1736 double const p_FR = -chi_S_L * p_cap_ip;
1737 // p_SR
1738 variables.solid_grain_pressure =
1739 p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
1740 auto const rho_SR =
1741 solid_phase.property(MPL::PropertyType::density)
1742 .template value<double>(variables, x_position, t, dt);
1743 *std::get<DrySolidDensity>(this->output_data_[ip]) = (1 - phi) * rho_SR;
1744
1745 {
1746 auto& SD = this->current_states_[ip];
1747 auto const& sigma_sw =
1748 std::get<ProcessLib::ThermoRichardsMechanics::
1749 ConstitutiveStress_StrainTemperature::
1750 SwellingDataStateful<DisplacementDim>>(SD)
1751 .sigma_sw;
1752 auto& eps_m =
1753 std::get<ProcessLib::ThermoRichardsMechanics::
1754 ConstitutiveStress_StrainTemperature::
1755 MechanicalStrainData<DisplacementDim>>(SD)
1756 .eps_m;
1757 eps_m.noalias() =
1758 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
1759 ? eps + C_el.inverse() * sigma_sw
1760 : eps;
1761 variables.mechanical_strain.emplace<
1762 MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
1763 eps_m);
1764 }
1765
1766 {
1767 auto& SD = this->current_states_[ip];
1768 auto const& SD_prev = this->prev_states_[ip];
1769 auto& sigma_eff =
1770 std::get<ProcessLib::ThermoRichardsMechanics::
1771 ConstitutiveStress_StrainTemperature::
1772 EffectiveStressData<DisplacementDim>>(SD);
1773 auto const& sigma_eff_prev = std::get<
1774 PrevState<ProcessLib::ThermoRichardsMechanics::
1775 ConstitutiveStress_StrainTemperature::
1776 EffectiveStressData<DisplacementDim>>>(
1777 SD_prev);
1778 auto const& eps_m =
1779 std::get<ProcessLib::ThermoRichardsMechanics::
1780 ConstitutiveStress_StrainTemperature::
1781 MechanicalStrainData<DisplacementDim>>(SD);
1782 auto const& eps_m_prev = std::get<
1783 PrevState<ProcessLib::ThermoRichardsMechanics::
1784 ConstitutiveStress_StrainTemperature::
1785 MechanicalStrainData<DisplacementDim>>>(
1786 SD_prev);
1787
1788 _ip_data[ip].updateConstitutiveRelation(
1789 variables, t, x_position, dt, temperature, sigma_eff,
1790 sigma_eff_prev, eps_m, eps_m_prev, this->solid_material_,
1791 this->material_states_[ip].material_state_variables);
1792 }
1793
1794 auto const& b = this->process_data_.specific_body_force;
1795
1796 // Compute the velocity
1797 auto const& dNdx_p = _ip_data[ip].dNdx_p;
1798 std::get<
1799 ProcessLib::ThermoRichardsMechanics::DarcyLawData<DisplacementDim>>(
1800 this->output_data_[ip])
1801 ->noalias() = -K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b;
1802
1803 saturation_avg += S_L;
1804 porosity_avg += phi;
1805 sigma_avg += sigma_eff;
1806 }
1807 saturation_avg /= n_integration_points;
1808 porosity_avg /= n_integration_points;
1809 sigma_avg /= n_integration_points;
1810
1811 (*this->process_data_.element_saturation)[this->element_.getID()] =
1812 saturation_avg;
1813 (*this->process_data_.element_porosity)[this->element_.getID()] =
1814 porosity_avg;
1815
1816 Eigen::Map<KV>(
1817 &(*this->process_data_.element_stresses)[this->element_.getID() *
1818 KV::RowsAtCompileTime]) =
1819 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_avg);
1820
1821 NumLib::interpolateToHigherOrderNodes<
1822 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1823 DisplacementDim>(this->element_, this->is_axially_symmetric_, p_L,
1824 *this->process_data_.pressure_interpolated);
1825}
1826} // namespace RichardsMechanics
1827} // namespace ProcessLib
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:26
MaterialPropertyLib::Medium::phase
Phase const & phase(std::size_t index) const
Definition Medium.cpp:33
MaterialPropertyLib::Medium::property
Property const & property(PropertyType const &p) const
Definition Medium.cpp:54
MaterialPropertyLib::Medium::hasProperty
bool hasProperty(PropertyType const &p) const
Definition Medium.cpp:70
MaterialPropertyLib::Phase::property
Property const & property(PropertyType const &p) const
Definition Phase.cpp:53
MaterialPropertyLib::Phase::hasProperty
bool hasProperty(PropertyType const &p) const
Definition Phase.cpp:69
MaterialPropertyLib::Property::value
virtual PropertyDataType value() const
Definition Property.cpp:76
ParameterLib::SpatialPosition::getElementID
std::optional< std::size_t > getElementID() const
Definition SpatialPosition.h:43
ParameterLib::SpatialPosition::setElementID
void setElementID(std::size_t element_id)
Definition SpatialPosition.h:59
ParameterLib::SpatialPosition::getIntegrationPoint
std::optional< unsigned > getIntegrationPoint() const
Definition SpatialPosition.h:44
ParameterLib::SpatialPosition::setIntegrationPoint
void setIntegrationPoint(unsigned integration_point)
Definition SpatialPosition.h:65
ProcessLib::BMatrixPolicyType::BMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition BMatrixPolicy.h:55
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::SymmetricTensor
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
Definition RichardsMechanicsFEM.h:76
ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler::GlobalDimMatrixType
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
Definition RichardsMechanicsFEM.h:60
MaterialPropertyLib::formEigenTensor< 3 >
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)
MathLib::KelvinVector::kelvinVectorToSymmetricTensor
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Definition KelvinVector.cpp:142
MathLib::KelvinVector::kelvin_vector_dimensions
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition KelvinVector.h:24
MathLib::KelvinVector::KelvinVectorType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition KelvinVector.h:47
MathLib::KelvinVector::KelvinMatrixType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
Definition KelvinVector.h:55
MathLib::createZeroedVector
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Definition EigenMapTools.h:153
MathLib::createZeroedMatrix
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:32
NumLib::detail::shapeFunctionInterpolate
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition Interpolation.h:27
NumLib::interpolateXCoordinate
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition InitShapeMatrices.h:71
NumLib::interpolateToHigherOrderNodes
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
Definition Interpolation.h:107
NumLib::initShapeMatrices
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)
Definition InitShapeMatrices.h:53
ProcessLib::LinearBMatrix::computeBMatrix
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.
Definition LinearBMatrix.h:46
ProcessLib::RichardsMechanics::StatefulData
std::tuple< StrainData< DisplacementDim >, ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature:: EffectiveStressData< DisplacementDim >, ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature:: SwellingDataStateful< DisplacementDim >, ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature:: MechanicalStrainData< DisplacementDim >, ProcessLib::ThermoRichardsMechanics::SaturationData, ProcessLib::ThermoRichardsMechanics::PorosityData, ProcessLib::ThermoRichardsMechanics::TransportPorosityData, MicroPressure, MicroSaturation > StatefulData
Data whose state must be tracked by the TRM process.
Definition ConstitutiveData.h:35
ProcessLib::RichardsMechanics::updateSwellingStressAndVolumetricStrain
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)
Definition RichardsMechanicsFEM-impl.h:35
ProcessLib::RichardsMechanics::StatefulDataPrev
ProcessLib::ConstitutiveRelations::PrevStateOf< StatefulData< DisplacementDim > > StatefulDataPrev
Definition ConstitutiveData.h:49
ProcessLib::RichardsMechanics::createConstitutiveModels
ConstitutiveModels< DisplacementDim > createConstitutiveModels(TRMProcessData const &process_data, MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
Definition ConstitutiveModels.h:24
ProcessLib::RichardsMechanics::ConstitutiveData
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.
Definition ConstitutiveData.h:60
ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature
Definition ConstitutiveData.h:34
EigenFixedShapeMatrixPolicy::NodalMatrixType
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
Definition ShapeMatrixPolicy.h:74
EigenFixedShapeMatrixPolicy::GlobalDimVectorType
VectorType< GlobalDim > GlobalDimVectorType
Definition ShapeMatrixPolicy.h:85
MaterialLib::Solids::MechanicsBase::getBulkModulus
virtual double getBulkModulus(double const, ParameterLib::SpatialPosition const &, KelvinMatrix const *const =nullptr) const
Definition MechanicsBase.h:136
ProcessLib::ConstitutiveRelations::PrevState
Represents a previous state of type T.
Definition Base.h:21
ProcessLib::RichardsMechanics::RichardsMechanicsProcessData::media_map
MaterialPropertyLib::MaterialSpatialDistributionMap media_map
Definition RichardsMechanicsProcessData.h:38