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