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