From e106fcf79bc5120d56b32e6388cdc9d5918484de Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Thu, 28 Mar 2024 17:02:27 +0100
Subject: [PATCH] [PL/TH2M] Extract porosity model
---
.../ConstitutiveModels.h | 7 ++
.../TH2M/ConstitutiveRelations/Porosity.cpp | 78 +++++++++++++++++++
.../TH2M/ConstitutiveRelations/Porosity.h | 28 ++++++-
ProcessLib/TH2M/TH2MFEM-impl.h | 41 ++--------
4 files changed, 120 insertions(+), 34 deletions(-)
create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index 7fe33f3683b..bb9847ebdc6 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -15,6 +15,7 @@
#include "MechanicalStrain.h"
#include "PermeabilityModel.h"
#include "PhaseTransitionModel.h"
+#include "Porosity.h"
#include "PureLiquidDensity.h"
#include "Saturation.h"
#include "SolidCompressibility.h"
@@ -58,6 +59,12 @@ struct ConstitutiveModels
PermeabilityModel<DisplacementDim> permeability_model;
PureLiquidDensityModel pure_liquid_density_model;
PhaseTransitionModel const& phase_transition_model;
+#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
+ PorosityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>
+#else // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
+ PorosityModel
+#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
+ porosity_model;
};
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp
new file mode 100644
index 00000000000..14533247d5e
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp
@@ -0,0 +1,78 @@
+/**
+ * \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
+ */
+
+#include "Porosity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+
+void PorosityModel::eval(SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ PorosityData& porosity_data,
+ PorosityDerivativeData& porosity_d_data) const
+{
+ MaterialPropertyLib::VariableArray variables;
+
+ auto const& mpl_porosity =
+ media_data.medium[MaterialPropertyLib::PropertyType::porosity];
+
+ porosity_data.phi =
+ mpl_porosity.template value<double>(variables, x_t.x, x_t.t, x_t.dt);
+
+ porosity_d_data.dphi_S_dT = -mpl_porosity.template dValue<double>(
+ variables, MaterialPropertyLib::Variable::temperature, x_t.x, x_t.t,
+ x_t.dt);
+}
+
+template <int DisplacementDim>
+void PorosityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>::eval(
+ SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ BiotData const& biot,
+ StrainData<DisplacementDim> const& strain_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ PorosityData& porosity_data,
+ PorosityDerivativeData& porosity_d_data) const
+{
+ MaterialPropertyLib::VariableArray variables;
+
+ auto const& mpl_porosity =
+ media_data.medium[MaterialPropertyLib::PropertyType::porosity];
+
+ double const phi_0 =
+ mpl_porosity.template value<double>(variables, x_t.x, x_t.t, x_t.dt);
+
+ static int const KelvinVectorSize =
+ MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
+ using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+ double const div_u = Invariants::trace(strain_data.eps);
+
+ double const phi_S =
+ (1. - phi_0) *
+ (1. + s_therm_exp_data.thermal_volume_strain - biot() * div_u);
+
+ porosity_data.phi = 1. - phi_S;
+
+ auto const dphi_0_dT = mpl_porosity.template dValue<double>(
+ variables, MaterialPropertyLib::Variable::temperature, x_t.x, x_t.t,
+ x_t.dt);
+
+ auto const dphi_S_0_dT = -dphi_0_dT;
+ porosity_d_data.dphi_S_dT =
+ dphi_S_0_dT *
+ (1. + s_therm_exp_data.thermal_volume_strain - biot() * div_u) +
+ (1. - phi_0) * s_therm_exp_data.beta_T_SR;
+}
+
+template struct PorosityModelNonConstantSolidPhaseVolumeFraction<2>;
+template struct PorosityModelNonConstantSolidPhaseVolumeFraction<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h
index f6b239a411d..5651fb31891 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h
@@ -10,7 +10,10 @@
#pragma once
#include "Base.h"
+#include "Biot.h"
+#include "ProcessLib/ConstitutiveRelations/StrainData.h"
#include "ProcessLib/Reflection/ReflectionData.h"
+#include "SolidThermalExpansion.h"
namespace ProcessLib::TH2M
{
@@ -19,7 +22,6 @@ namespace ConstitutiveRelations
struct PorosityDerivativeData
{
double dphi_S_dT = nan;
- double phi_0 = nan; // Initial porosity.
};
struct PorosityData
@@ -35,5 +37,29 @@ struct PorosityData
}
};
+struct PorosityModel
+{
+ void eval(SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ PorosityData& porosity_data,
+ PorosityDerivativeData& porosity_d_data) const;
+};
+
+template <int DisplacementDim>
+struct PorosityModelNonConstantSolidPhaseVolumeFraction
+{
+ void eval(
+ SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ BiotData const& biot,
+ StrainData<DisplacementDim> const& strain_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ PorosityData& porosity_data,
+ PorosityDerivativeData& porosity_d_data) const;
+};
+
+extern template struct PorosityModelNonConstantSolidPhaseVolumeFraction<2>;
+extern template struct PorosityModelNonConstantSolidPhaseVolumeFraction<3>;
+
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index e6192dde969..93bb5ae7f61 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -232,6 +232,13 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ip_out.vapour_pressure_data, current_state.constituent_density_data,
ip_cv.phase_transition_data);
+ models.porosity_model.eval({pos, t, dt}, media_data,
+#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
+ ip_cv.biot_data, ip_out.eps_data,
+ ip_cv.s_therm_exp_data,
+#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
+ ip_out.porosity_data, ip_cv.porosity_d_data);
+
MPL::VariableArray vars;
MPL::VariableArray vars_prev;
vars.temperature = T;
@@ -245,23 +252,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
vars.liquid_saturation = current_state.S_L_data.S_L;
vars_prev.liquid_saturation = prev_state.S_L_data->S_L;
- ip_cv.porosity_d_data.phi_0 =
- medium.property(MPL::PropertyType::porosity)
- .template value<double>(vars, pos, t, dt);
-
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- auto const& m = Invariants::identity2;
- double const div_u = m.transpose() * eps;
-
- double const phi_S =
- (1. - ip_cv.porosity_d_data.phi_0) *
- (1. + ip_data.thermal_volume_strain - ip_cv.biot_data() * div_u);
-#else // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- double const phi_S = (1. - ip_cv.porosity_d_data.phi_0);
-#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
-
- // porosity
- ip_out.porosity_data.phi = 1. - phi_S;
vars.porosity = ip_out.porosity_data.phi;
// solid phase density
@@ -284,6 +274,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
current_state.S_L_data.S_L * ip_out.porosity_data.phi;
auto const phi_G =
(1. - current_state.S_L_data.S_L) * ip_out.porosity_data.phi;
+ double const phi_S = 1. - ip_out.porosity_data.phi;
// thermal conductivity
ip_data.lambda = MaterialPropertyLib::formEigenTensor<DisplacementDim>(
@@ -369,22 +360,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
#endif
;
- // porosity
- auto const dphi_0_dT =
- medium[MPL::PropertyType::porosity].template dValue<double>(
- vars, MPL::Variable::temperature, pos, t, dt);
-
- auto const dphi_S_0_dT = -dphi_0_dT;
- ip_cv.porosity_d_data.dphi_S_dT =
- dphi_S_0_dT
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- * (1. + ip_cv.s_therm_exp_data.thermal_volume_strain -
- ip_cv.biot_data() * div_u) +
- (1. - ip_cv.porosity_d_data.phi_0) *
- ip_cv.s_therm_exp_data.beta_T_SR
-#endif
- ;
-
ip_cv.drho_u_eff_dT =
phi_G * c.drho_GR_dT * c.uG +
phi_G * ip_out.fluid_density_data.rho_GR * c.du_G_dT +
--
GitLab