From 2fe391a69a66b78306a0738af01e372276faa0d6 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Tue, 28 May 2024 12:11:56 +0200
Subject: [PATCH] [PL/RM] Extracted some constitutive setting evaluations
... from assembleWithJacobian() into a separate method.
I added some constitutive data structs and reused some from TRM.
---
.../ConstitutiveRelations/Base.h | 9 +-
.../ConstitutiveRelations/ConstitutiveData.h | 23 +-
.../ConstitutiveSetting.cpp | 24 +-
.../ConstitutiveSetting.h | 3 +-
.../ConstitutiveRelations/Density.h | 19 +
.../EffectivePorePressure.h | 19 +
.../ConstitutiveRelations/LiquidDensity.h | 18 +
.../SaturationSecantDerivative.h | 21 +
.../ConstitutiveRelations/StiffnessTensor.h | 20 +
.../RichardsMechanicsFEM-impl.h | 424 +++++++++++-------
.../RichardsMechanics/RichardsMechanicsFEM.h | 10 +
11 files changed, 410 insertions(+), 180 deletions(-)
create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/Density.h
create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/EffectivePorePressure.h
create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/LiquidDensity.h
create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/SaturationSecantDerivative.h
create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/StiffnessTensor.h
diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/Base.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/Base.h
index be08c171856..a70b6227fc9 100644
--- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/Base.h
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/Base.h
@@ -10,6 +10,7 @@
#pragma once
+#include "BaseLib/StrongType.h"
#include "MaterialLib/MPL/Medium.h"
#include "MathLib/KelvinVector.h"
#include "ParameterLib/SpatialPosition.h"
@@ -63,13 +64,7 @@ struct MediaData
MaterialPropertyLib::Phase const& solid;
};
-template <int DisplacementDim>
-struct TemperatureData
-{
- double T;
- double T_prev;
- Eigen::Vector<double, DisplacementDim> grad_T;
-};
+using TemperatureData = BaseLib::StrongType<double, struct TemperatureDataTag>;
template <int DisplacementDim>
struct CapillaryPressureData
diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h
index 8b750271bdb..e909848e70b 100644
--- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h
@@ -10,7 +10,17 @@
#pragma once
+#include "Density.h"
+#include "LiquidDensity.h"
#include "ProcessLib/ConstitutiveRelations/Base.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Biot.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/LiquidViscosity.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityData.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Porosity.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Saturation.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/SolidCompressibilityData.h"
+#include "SaturationSecantDerivative.h"
+#include "StiffnessTensor.h"
namespace ProcessLib::RichardsMechanics
{
@@ -29,7 +39,18 @@ using OutputData = std::tuple<>;
/// Data that is needed for the equation system assembly.
template <int DisplacementDim>
-using ConstitutiveData = std::tuple<>;
+using ConstitutiveData = std::tuple<
+ // TODO (CL) check if all that data should stay here
+ 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>;
/// Data that stores intermediate values, which are not needed outside the
/// constitutive setting.
diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.cpp b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.cpp
index 836ae0d3b7f..63dd576bd85 100644
--- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.cpp
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.cpp
@@ -14,19 +14,19 @@ namespace ProcessLib::RichardsMechanics
{
template <int DisplacementDim>
void ConstitutiveSetting<DisplacementDim>::eval(
- ConstitutiveModels<DisplacementDim>& models, double const t,
- double const dt, ParameterLib::SpatialPosition const& x_position,
- MaterialPropertyLib::Medium const& medium,
- TemperatureData<DisplacementDim> const& T_data,
- CapillaryPressureData<DisplacementDim> const& p_cap_data,
- KelvinVector<DisplacementDim> const& eps_arg,
- StatefulData<DisplacementDim>& state,
- StatefulDataPrev<DisplacementDim> const& prev_state,
+ ConstitutiveModels<DisplacementDim>& /*models*/, double const /*t*/,
+ double const /*dt*/, ParameterLib::SpatialPosition const& /*x_position*/,
+ MaterialPropertyLib::Medium const& /*medium*/,
+ TemperatureData const /*T_data*/,
+ CapillaryPressureData<DisplacementDim> const& /*p_cap_data*/,
+ KelvinVector<DisplacementDim> const& /*eps_arg*/,
+ StatefulData<DisplacementDim>& /*state*/,
+ StatefulDataPrev<DisplacementDim> const& /*prev_state*/,
ProcessLib::ThermoRichardsMechanics::MaterialStateData<DisplacementDim>&
- mat_state,
- ConstitutiveTempData<DisplacementDim>& tmp,
- OutputData<DisplacementDim>& out,
- ConstitutiveData<DisplacementDim>& cd) const
+ /*mat_state*/,
+ ConstitutiveTempData<DisplacementDim>& /*tmp*/,
+ OutputData<DisplacementDim>& /*out*/,
+ ConstitutiveData<DisplacementDim>& /*cd*/) const
{
}
diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.h
index 3dbae0f6195..8e10f44ec91 100644
--- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.h
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.h
@@ -24,8 +24,7 @@ struct ConstitutiveSetting
void eval(
ConstitutiveModels<DisplacementDim>& models, double const t,
double const dt, ParameterLib::SpatialPosition const& x_position,
- MaterialPropertyLib::Medium const& medium,
- TemperatureData<DisplacementDim> const& T_data,
+ MaterialPropertyLib::Medium const& medium, TemperatureData const T_data,
CapillaryPressureData<DisplacementDim> const& p_cap_data,
KelvinVector<DisplacementDim> const& eps_arg,
StatefulData<DisplacementDim>& state,
diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/Density.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/Density.h
new file mode 100644
index 00000000000..f53cd357b24
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/Density.h
@@ -0,0 +1,19 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "BaseLib/StrongType.h"
+
+namespace ProcessLib::RichardsMechanics
+{
+// effective density of both the solid and fluid phases
+using Density = BaseLib::StrongType<double, struct DensityTag>;
+} // namespace ProcessLib::RichardsMechanics
diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/EffectivePorePressure.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/EffectivePorePressure.h
new file mode 100644
index 00000000000..5df96875c01
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/EffectivePorePressure.h
@@ -0,0 +1,19 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "BaseLib/StrongType.h"
+
+namespace ProcessLib::RichardsMechanics
+{
+using EffectivePorePressure =
+ BaseLib::StrongType<double, struct EffectivePorePressureTag>;
+} // namespace ProcessLib::RichardsMechanics
diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/LiquidDensity.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/LiquidDensity.h
new file mode 100644
index 00000000000..05b11b32c31
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/LiquidDensity.h
@@ -0,0 +1,18 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "BaseLib/StrongType.h"
+
+namespace ProcessLib::RichardsMechanics
+{
+using LiquidDensity = BaseLib::StrongType<double, struct LiquidDensityTag>;
+} // namespace ProcessLib::RichardsMechanics
diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/SaturationSecantDerivative.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/SaturationSecantDerivative.h
new file mode 100644
index 00000000000..0e70b4cbc0f
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/SaturationSecantDerivative.h
@@ -0,0 +1,21 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "Base.h"
+
+namespace ProcessLib::RichardsMechanics
+{
+struct SaturationSecantDerivative
+{
+ double DeltaS_L_Deltap_cap = nan;
+};
+} // namespace ProcessLib::RichardsMechanics
diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/StiffnessTensor.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/StiffnessTensor.h
new file mode 100644
index 00000000000..46b272c2d95
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/StiffnessTensor.h
@@ -0,0 +1,20 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "Base.h"
+
+namespace ProcessLib::RichardsMechanics
+{
+template <int DisplacementDim>
+using StiffnessTensor = BaseLib::StrongType<KelvinMatrix<DisplacementDim>,
+ struct StiffnessTensorTag>;
+} // namespace ProcessLib::RichardsMechanics
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 87647190529..54a0d49d9be 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -685,6 +685,216 @@ void RichardsMechanicsLocalAssembler<
}
}
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ int DisplacementDim>
+void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure, DisplacementDim>::
+ assembleWithJacobianEvalConstitutiveSetting(
+ double const t, double const dt,
+ ParameterLib::SpatialPosition const& x_position,
+ RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure,
+ DisplacementDim>::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)
+{
+ auto const& liquid_phase = medium->phase("AqueousLiquid");
+ auto const& solid_phase = medium->phase("Solid");
+
+ auto const& identity2 = MathLib::KelvinVector::Invariants<
+ MathLib::KelvinVector::kelvin_vector_dimensions(
+ DisplacementDim)>::identity2;
+
+ double const temperature = T_data();
+ double const p_cap_ip = p_cap_data.p_cap;
+ double const p_cap_prev_ip = p_cap_data.p_cap_prev;
+
+ auto& eps = ip_data.eps;
+ auto& eps_m = ip_data.eps_m;
+ auto& S_L = ip_data.saturation;
+ auto const S_L_prev = ip_data.saturation_prev;
+ auto const alpha =
+ medium->property(MPL::PropertyType::biot_coefficient)
+ .template value<double>(variables, x_position, t, dt);
+ *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD) = alpha;
+
+ auto const C_el =
+ ip_data.computeElasticTangentStiffness(t, x_position, dt, temperature);
+
+ auto const beta_SR = (1 - alpha) / ip_data.solid_material.getBulkModulus(
+ t, x_position, &C_el);
+ variables.grain_compressibility = beta_SR;
+ std::get<ProcessLib::ThermoRichardsMechanics::SolidCompressibilityData>(CD)
+ .beta_SR = beta_SR;
+
+ auto const rho_LR =
+ liquid_phase.property(MPL::PropertyType::density)
+ .template value<double>(variables, x_position, t, dt);
+ variables.density = rho_LR;
+ *std::get<LiquidDensity>(CD) = rho_LR;
+
+ S_L = medium->property(MPL::PropertyType::saturation)
+ .template value<double>(variables, x_position, t, dt);
+ variables.liquid_saturation = S_L;
+ variables_prev.liquid_saturation = S_L_prev;
+
+ // tangent derivative for Jacobian
+ double const dS_L_dp_cap =
+ medium->property(MPL::PropertyType::saturation)
+ .template dValue<double>(variables,
+ MPL::Variable::capillary_pressure,
+ x_position, t, dt);
+ std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(CD)
+ .dS_L_dp_cap = dS_L_dp_cap;
+ // secant derivative from time discretization for storage
+ // use tangent, if secant is not available
+ double const DeltaS_L_Deltap_cap =
+ (p_cap_ip == p_cap_prev_ip)
+ ? dS_L_dp_cap
+ : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
+ std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap =
+ DeltaS_L_Deltap_cap;
+
+ auto const chi = [medium, x_position, t, dt](double const S_L)
+ {
+ MPL::VariableArray vs;
+ vs.liquid_saturation = S_L;
+ return medium->property(MPL::PropertyType::bishops_effective_stress)
+ .template value<double>(vs, x_position, t, dt);
+ };
+ double const chi_S_L = chi(S_L);
+ std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD).chi_S_L =
+ chi_S_L;
+ double const chi_S_L_prev = chi(S_L_prev);
+ std::get<PrevState<ProcessLib::ThermoRichardsMechanics::BishopsData>>(CD)
+ ->chi_S_L = chi_S_L_prev;
+
+ auto const dchi_dS_L =
+ medium->property(MPL::PropertyType::bishops_effective_stress)
+ .template dValue<double>(
+ variables, MPL::Variable::liquid_saturation, x_position, t, dt);
+ std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD).dchi_dS_L =
+ dchi_dS_L;
+
+ double const p_FR = -chi_S_L * p_cap_ip;
+ variables.effective_pore_pressure = p_FR;
+ variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
+
+ // Set volumetric strain rate for the general case without swelling.
+ variables.volumetric_strain = Invariants::trace(ip_data.eps);
+ // TODO (CL) changed that, using eps_prev for the moment, not B * u_prev
+ // variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
+ variables_prev.volumetric_strain = Invariants::trace(ip_data.eps_prev);
+
+ auto& phi = ip_data.porosity;
+ { // Porosity update
+
+ variables_prev.porosity = ip_data.porosity_prev;
+ phi = medium->property(MPL::PropertyType::porosity)
+ .template value<double>(variables, variables_prev, x_position,
+ t, dt);
+ variables.porosity = phi;
+ }
+ std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi = phi;
+
+ if (alpha < phi)
+ {
+ OGS_FATAL(
+ "RichardsMechanics: Biot-coefficient {} is smaller than "
+ "porosity {} in element/integration point {}/{}.",
+ alpha, phi, _element.getID(), -1 /*ip*/ /* TODO (CL) re-enable */);
+ }
+
+ auto const mu = liquid_phase.property(MPL::PropertyType::viscosity)
+ .template value<double>(variables, x_position, t, dt);
+ std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(CD)
+ .viscosity = mu;
+
+ // Swelling and possibly volumetric strain rate update.
+ updateSwellingStressAndVolumetricStrain<DisplacementDim>(
+ ip_data, *medium, solid_phase, C_el, rho_LR, mu,
+ _process_data.micro_porosity_parameters, alpha, phi, p_cap_ip,
+ variables, variables_prev, x_position, t, dt);
+
+ if (medium->hasProperty(MPL::PropertyType::transport_porosity))
+ {
+ if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
+ {
+ variables_prev.transport_porosity = ip_data.transport_porosity_prev;
+
+ ip_data.transport_porosity =
+ medium->property(MPL::PropertyType::transport_porosity)
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
+ variables.transport_porosity = ip_data.transport_porosity;
+ }
+ }
+ else
+ {
+ variables.transport_porosity = phi;
+ }
+
+ // Set mechanical variables for the intrinsic permeability model
+ // For stress dependent permeability.
+ {
+ auto const sigma_total =
+ (ip_data.sigma_eff + alpha * p_FR * identity2).eval();
+ // For stress dependent permeability.
+ variables.total_stress.emplace<SymmetricTensor>(
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
+ }
+
+ variables.equivalent_plastic_strain =
+ ip_data.material_state_variables->getEquivalentPlasticStrain();
+
+ double const k_rel =
+ medium->property(MPL::PropertyType::relative_permeability)
+ .template value<double>(variables, x_position, t, dt);
+
+ auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
+ medium->property(MPL::PropertyType::permeability)
+ .value(variables, x_position, t, dt));
+
+ std::get<
+ ProcessLib::ThermoRichardsMechanics::PermeabilityData<DisplacementDim>>(
+ CD)
+ .k_rel = k_rel;
+ std::get<
+ ProcessLib::ThermoRichardsMechanics::PermeabilityData<DisplacementDim>>(
+ CD)
+ .Ki = K_intrinsic;
+
+ //
+ // displacement equation, displacement part
+ //
+ auto& sigma_sw = ip_data.sigma_sw;
+
+ eps_m.noalias() =
+ solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
+ ? eps + C_el.inverse() * sigma_sw
+ : eps;
+ variables.mechanical_strain
+ .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
+ eps_m);
+
+ auto C = ip_data.updateConstitutiveRelation(variables, t, x_position, dt,
+ temperature);
+
+ *std::get<StiffnessTensor<DisplacementDim>>(CD) = std::move(C);
+
+ // p_SR
+ variables.solid_grain_pressure =
+ p_FR - ip_data.sigma_eff.dot(identity2) / (3 * (1 - phi));
+ auto const rho_SR =
+ solid_phase.property(MPL::PropertyType::density)
+ .template value<double>(variables, x_position, t, dt);
+
+ double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
+ *std::get<Density>(CD) = rho;
+}
+
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
int DisplacementDim>
void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
@@ -775,7 +985,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
_integration_method.getNumberOfPoints();
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
- [[maybe_unused]] ConstitutiveData<DisplacementDim> CD;
+ ConstitutiveData<DisplacementDim> CD;
[[maybe_unused]] auto models = createConstitutiveModels(
_process_data, _ip_data[ip].solid_material);
@@ -816,166 +1026,25 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables.temperature = temperature;
auto& eps = _ip_data[ip].eps;
- auto& eps_m = _ip_data[ip].eps_m;
eps.noalias() = B * u;
- auto const& sigma_eff = _ip_data[ip].sigma_eff;
- auto& S_L = _ip_data[ip].saturation;
- auto const S_L_prev = _ip_data[ip].saturation_prev;
- auto const alpha =
- medium->property(MPL::PropertyType::biot_coefficient)
- .template value<double>(variables, x_position, t, dt);
- auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
- t, x_position, dt, temperature);
-
- auto const beta_SR =
- (1 - alpha) /
- _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
- variables.grain_compressibility = beta_SR;
-
- auto const rho_LR =
- liquid_phase.property(MPL::PropertyType::density)
- .template value<double>(variables, x_position, t, dt);
- variables.density = rho_LR;
-
- auto const& b = _process_data.specific_body_force;
-
- S_L = medium->property(MPL::PropertyType::saturation)
- .template value<double>(variables, x_position, t, dt);
- variables.liquid_saturation = S_L;
- variables_prev.liquid_saturation = S_L_prev;
-
- // tangent derivative for Jacobian
- double const dS_L_dp_cap =
- medium->property(MPL::PropertyType::saturation)
- .template dValue<double>(variables,
- MPL::Variable::capillary_pressure,
- x_position, t, dt);
- // secant derivative from time discretization for storage
- // use tangent, if secant is not available
- double const DeltaS_L_Deltap_cap =
- (p_cap_ip == p_cap_prev_ip)
- ? dS_L_dp_cap
- : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
-
- auto const chi = [medium, x_position, t, dt](double const S_L)
- {
- MPL::VariableArray vs;
- vs.liquid_saturation = S_L;
- return medium->property(MPL::PropertyType::bishops_effective_stress)
- .template value<double>(vs, x_position, t, dt);
- };
- double const chi_S_L = chi(S_L);
- double const chi_S_L_prev = chi(S_L_prev);
-
- double const p_FR = -chi_S_L * p_cap_ip;
- variables.effective_pore_pressure = p_FR;
- variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
-
- // Set volumetric strain rate for the general case without swelling.
- variables.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
- variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
-
- auto& phi = _ip_data[ip].porosity;
- { // Porosity update
-
- variables_prev.porosity = _ip_data[ip].porosity_prev;
- phi = medium->property(MPL::PropertyType::porosity)
- .template value<double>(variables, variables_prev,
- x_position, t, dt);
- variables.porosity = phi;
- }
-
- if (alpha < phi)
- {
- OGS_FATAL(
- "RichardsMechanics: Biot-coefficient {} is smaller than "
- "porosity {} in element/integration point {}/{}.",
- alpha, phi, _element.getID(), ip);
- }
-
- auto const mu =
- liquid_phase.property(MPL::PropertyType::viscosity)
- .template value<double>(variables, x_position, t, dt);
-
- // Swelling and possibly volumetric strain rate update.
- updateSwellingStressAndVolumetricStrain<DisplacementDim>(
- _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu,
- _process_data.micro_porosity_parameters, alpha, phi, p_cap_ip,
- variables, variables_prev, x_position, t, dt);
-
- if (medium->hasProperty(MPL::PropertyType::transport_porosity))
- {
- if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
- {
- variables_prev.transport_porosity =
- _ip_data[ip].transport_porosity_prev;
-
- _ip_data[ip].transport_porosity =
- medium->property(MPL::PropertyType::transport_porosity)
- .template value<double>(variables, variables_prev,
- x_position, t, dt);
- variables.transport_porosity = _ip_data[ip].transport_porosity;
- }
- }
- else
- {
- variables.transport_porosity = phi;
- }
-
- double const k_rel =
- medium->property(MPL::PropertyType::relative_permeability)
- .template value<double>(variables, x_position, t, dt);
-
- // Set mechanical variables for the intrinsic permeability model
- // For stress dependent permeability.
- {
- auto const sigma_total =
- (_ip_data[ip].sigma_eff + alpha * p_FR * identity2).eval();
- // For stress dependent permeability.
- variables.total_stress.emplace<SymmetricTensor>(
- MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
- sigma_total));
- }
-
- variables.equivalent_plastic_strain =
- _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
-
- auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
- medium->property(MPL::PropertyType::permeability)
- .value(variables, x_position, t, dt));
-
- GlobalDimMatrixType const rho_Ki_over_mu = K_intrinsic * rho_LR / mu;
-
- //
- // displacement equation, displacement part
- //
- auto& sigma_sw = _ip_data[ip].sigma_sw;
-
- eps_m.noalias() =
- solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
- ? eps + C_el.inverse() * sigma_sw
- : eps;
- variables.mechanical_strain
- .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
- eps_m);
-
- auto C = _ip_data[ip].updateConstitutiveRelation(
- variables, t, x_position, dt, temperature);
+ assembleWithJacobianEvalConstitutiveSetting(
+ t, dt, x_position, _ip_data[ip], variables, variables_prev, medium,
+ TemperatureData{temperature},
+ CapillaryPressureData<DisplacementDim>{
+ p_cap_ip, p_cap_prev_ip,
+ Eigen::Vector<double, DisplacementDim>::Zero()},
+ CD);
+ auto const& C = *std::get<StiffnessTensor<DisplacementDim>>(CD);
local_Jac
.template block<displacement_size, displacement_size>(
displacement_index, displacement_index)
.noalias() += B.transpose() * C * B * w;
- // p_SR
- variables.solid_grain_pressure =
- p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
- auto const rho_SR =
- solid_phase.property(MPL::PropertyType::density)
- .template value<double>(variables, x_position, t, dt);
-
- double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
+ auto const& sigma_eff = _ip_data[ip].sigma_eff;
+ double const rho = *std::get<Density>(CD);
+ auto const& b = _process_data.specific_body_force;
local_rhs.template segment<displacement_size>(displacement_index)
.noalias() -=
(B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
@@ -983,13 +1052,21 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
//
// displacement equation, pressure part
//
+
+ double const alpha =
+ *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD);
+ double const chi_S_L =
+ std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
+ .chi_S_L;
Kup.noalias() += B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
- auto const dchi_dS_L =
- medium->property(MPL::PropertyType::bishops_effective_stress)
- .template dValue<double>(variables,
- MPL::Variable::liquid_saturation,
- x_position, t, dt);
+ double const dchi_dS_L =
+ std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
+ .dchi_dS_L;
+ double const dS_L_dp_cap =
+ std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(
+ CD)
+ .dS_L_dp_cap;
local_Jac
.template block<displacement_size, pressure_size>(
displacement_index, pressure_index)
@@ -997,6 +1074,9 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
(chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
identity2 * N_p * w;
+ double const phi =
+ std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi;
+ double const rho_LR = *std::get<LiquidDensity>(CD);
local_Jac
.template block<displacement_size, pressure_size>(
displacement_index, pressure_index)
@@ -1033,8 +1113,12 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
//
// pressure equation, displacement part.
//
+ double const S_L = _ip_data[ip].saturation;
if (_process_data.explicit_hm_coupling_in_unsaturated_zone)
{
+ double const chi_S_L_prev = std::get<PrevState<
+ ProcessLib::ThermoRichardsMechanics::BishopsData>>(CD)
+ ->chi_S_L;
Kpu.noalias() += N_p.transpose() * chi_S_L_prev * rho_LR * alpha *
identity2.transpose() * B * w;
}
@@ -1047,6 +1131,22 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
//
// pressure equation, pressure part.
//
+
+ double const k_rel =
+ std::get<ProcessLib::ThermoRichardsMechanics::PermeabilityData<
+ DisplacementDim>>(CD)
+ .k_rel;
+ auto const& K_intrinsic =
+ std::get<ProcessLib::ThermoRichardsMechanics::PermeabilityData<
+ DisplacementDim>>(CD)
+ .Ki;
+ double const mu =
+ std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(
+ CD)
+ .viscosity;
+
+ GlobalDimMatrixType const rho_Ki_over_mu = K_intrinsic * rho_LR / mu;
+
laplace_p.noalias() +=
dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
@@ -1057,6 +1157,11 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
MPL::Variable::liquid_phase_pressure,
x_position, t, dt);
+ double const beta_SR =
+ std::get<
+ ProcessLib::ThermoRichardsMechanics::SolidCompressibilityData>(
+ CD)
+ .beta_SR;
double const a0 = (alpha - phi) * beta_SR;
double const specific_storage_a_p = S_L * (phi * beta_LR + S_L * a0);
double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
@@ -1069,6 +1174,8 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
storage_p_a_p.noalias() +=
N_p.transpose() * rho_LR * specific_storage_a_p * N_p * w;
+ double const DeltaS_L_Deltap_cap =
+ std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap;
storage_p_a_S.noalias() -= N_p.transpose() * rho_LR *
specific_storage_a_S * DeltaS_L_Deltap_cap *
N_p * w;
@@ -1079,6 +1186,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
.noalias() += N_p.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
rho_LR * dspecific_storage_a_p_dp_cap * N_p * w;
+ double const S_L_prev = _ip_data[ip].saturation_prev;
storage_p_a_S_Jpp.noalias() -=
N_p.transpose() * rho_LR *
((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
index 7722a8a6cb2..1ad821a245e 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
@@ -13,6 +13,8 @@
#include <memory>
#include <vector>
+#include "ConstitutiveRelations/Base.h"
+#include "ConstitutiveRelations/ConstitutiveData.h"
#include "IntegrationPointData.h"
#include "LocalAssemblerInterface.h"
#include "MaterialLib/MPL/VariableType.h"
@@ -328,6 +330,14 @@ private:
getMaterialStateVariablesAt(unsigned integration_point) const override;
private:
+ 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);
+
RichardsMechanicsProcessData<DisplacementDim>& _process_data;
std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
--
GitLab