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