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