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