diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Advection.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Advection.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a236f888782dd03735b88c69f5727c6629682229
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Advection.cpp
@@ -0,0 +1,75 @@
+/**
+ * \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 "Advection.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void AdvectionModel<DisplacementDim>::eval(
+ ConstituentDensityData const& constituent_density_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ PureLiquidDensityData const& rho_W_LR,
+ ViscosityData const& viscosity_data,
+ AdvectionData<DisplacementDim>& advection_data) const
+{
+ GlobalDimMatrix<DisplacementDim> const k_over_mu_G =
+ permeability_data.Ki * permeability_data.k_rel_G / viscosity_data.mu_GR;
+ GlobalDimMatrix<DisplacementDim> const k_over_mu_L =
+ permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR;
+
+ advection_data.advection_C_G =
+ constituent_density_data.rho_C_GR * k_over_mu_G;
+ advection_data.advection_C_L =
+ constituent_density_data.rho_C_LR * k_over_mu_L;
+ advection_data.advection_W_G =
+ constituent_density_data.rho_W_GR * k_over_mu_G;
+ advection_data.advection_W_L = rho_W_LR() * k_over_mu_L;
+}
+
+template <int DisplacementDim>
+void AdvectionModel<DisplacementDim>::dEval(
+ ConstituentDensityData const& constituent_density_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ ViscosityData const& viscosity_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ PhaseTransitionData const& phase_transition_data,
+ AdvectionDerivativeData<DisplacementDim>& advection_d_data) const
+{
+ GlobalDimMatrix<DisplacementDim> const k_over_mu_G =
+ permeability_data.Ki * permeability_data.k_rel_G / viscosity_data.mu_GR;
+ GlobalDimMatrix<DisplacementDim> const k_over_mu_L =
+ permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR;
+
+ GlobalDimMatrix<DisplacementDim> const dk_over_mu_G_dp_cap =
+ permeability_data.Ki * permeability_data.dk_rel_G_dS_L * dS_L_dp_cap() /
+ viscosity_data.mu_GR;
+ GlobalDimMatrix<DisplacementDim> const dk_over_mu_L_dp_cap =
+ permeability_data.Ki * permeability_data.dk_rel_L_dS_L * dS_L_dp_cap() /
+ viscosity_data.mu_LR;
+
+ advection_d_data.dadvection_C_dp_GR =
+ phase_transition_data.drho_C_GR_dp_GR * k_over_mu_G
+ // + rhoCGR * (dk_over_mu_G_dp_GR = 0)
+ // + rhoCLR * (dk_over_mu_L_dp_GR = 0)
+ + phase_transition_data.drho_C_LR_dp_GR * k_over_mu_L;
+
+ advection_d_data.dadvection_C_dp_cap =
+ //(drho_C_GR_dp_cap = 0) * k_over_mu_G
+ constituent_density_data.rho_C_GR * dk_over_mu_G_dp_cap +
+ (-phase_transition_data.drho_C_LR_dp_LR) * k_over_mu_L +
+ constituent_density_data.rho_C_LR * dk_over_mu_L_dp_cap;
+}
+
+template struct AdvectionModel<2>;
+template struct AdvectionModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Advection.h b/ProcessLib/TH2M/ConstitutiveRelations/Advection.h
new file mode 100644
index 0000000000000000000000000000000000000000..316f52b82abc11c4cfb4d8932fe7312ffcc4278f
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Advection.h
@@ -0,0 +1,61 @@
+/**
+ * \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"
+#include "ConstitutiveDensity.h"
+#include "PermeabilityData.h"
+#include "PhaseTransitionData.h"
+#include "PureLiquidDensity.h"
+#include "Saturation.h"
+#include "Viscosity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+struct AdvectionData
+{
+ GlobalDimMatrix<DisplacementDim> advection_C_G;
+ GlobalDimMatrix<DisplacementDim> advection_C_L;
+ GlobalDimMatrix<DisplacementDim> advection_W_G;
+ GlobalDimMatrix<DisplacementDim> advection_W_L;
+};
+
+template <int DisplacementDim>
+struct AdvectionDerivativeData
+{
+ GlobalDimMatrix<DisplacementDim> dadvection_C_dp_GR;
+ GlobalDimMatrix<DisplacementDim> dadvection_C_dp_cap;
+};
+
+template <int DisplacementDim>
+struct AdvectionModel
+{
+ void eval(ConstituentDensityData const& constituent_density_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ PureLiquidDensityData const& rho_W_LR,
+ ViscosityData const& viscosity_data,
+ AdvectionData<DisplacementDim>& advection_data) const;
+
+ void dEval(
+ ConstituentDensityData const& constituent_density_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ ViscosityData const& viscosity_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ PhaseTransitionData const& phase_transition_data,
+ AdvectionDerivativeData<DisplacementDim>& advection_d_data) const;
+};
+
+extern template struct AdvectionModel<2>;
+extern template struct AdvectionModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Base.h b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
index 0667e3730d94f05cdd4c2b13063fe78061f6a1e1..70c0467e6cf95b44b5e6bd4da68ca8f3b574574e 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Base.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
@@ -59,6 +59,23 @@ using ReferenceTemperatureData =
using GasPressureData = BaseLib::StrongType<double, struct GasPressureTag>;
using CapillaryPressureData =
BaseLib::StrongType<double, struct CapillaryPressureTag>;
+template <int DisplacementDim>
+using GasPressureGradientData =
+ BaseLib::StrongType<GlobalDimVector<DisplacementDim>,
+ struct GasPressureGradientTag>;
+template <int DisplacementDim>
+using CapillaryPressureGradientData =
+ BaseLib::StrongType<GlobalDimVector<DisplacementDim>,
+ struct CapillaryPressureGradientTag>;
+template <int DisplacementDim>
+using TemperatureGradientData =
+ BaseLib::StrongType<GlobalDimVector<DisplacementDim>,
+ struct TemperatureGradientTag>;
+
+template <int DisplacementDim>
+using SpecificBodyForceData =
+ BaseLib::StrongType<GlobalDimVector<DisplacementDim>,
+ struct SpecificBodyForceTag>;
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/CEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/CEquation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6242c36efd48b3f4121e8d24e4ca349efd556223
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/CEquation.cpp
@@ -0,0 +1,467 @@
+/**
+ * \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 "CEquation.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void FC1Model<DisplacementDim>::eval(
+ AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ FC1Data<DisplacementDim>& fC_1) const
+{
+ fC_1.A = advection_data.advection_C_G * fluid_density_data.rho_GR +
+ advection_data.advection_C_L * fluid_density_data.rho_LR;
+}
+
+template struct FC1Model<2>;
+template struct FC1Model<3>;
+
+void FC2aModel::eval(BiotData const biot_data,
+ CapillaryPressureData const pCap,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const beta_p_SR,
+ FC2aData& fC_2a) const
+{
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ double const rho_C_FR = S_G * constituent_density_data.rho_C_GR +
+ S_L * constituent_density_data.rho_C_LR;
+ fC_2a.a =
+ porosity_data.phi * (constituent_density_data.rho_C_LR -
+ constituent_density_data.rho_C_GR) -
+ rho_C_FR * pCap() * (biot_data() - porosity_data.phi) * beta_p_SR();
+}
+
+void FC2aModel::dEval(BiotData const& biot_data,
+ CapillaryPressureData const pCap,
+ ConstituentDensityData const& constituent_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PorosityDerivativeData const& porosity_d_data,
+ SaturationData const& S_L_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ SolidCompressibilityData const& beta_p_SR,
+ FC2aDerivativeData& dfC_2a) const
+{
+ double const S_L = S_L_data.S_L;
+ double const S_G = 1. - S_L;
+
+ double const drho_C_FR_dp_GR =
+ /*(dS_G_dp_GR = 0) * constituent_density_data.rho_C_GR +*/
+ S_G * phase_transition_data.drho_C_GR_dp_GR +
+ /*(dS_L_dp_GR = 0) * constituent_density_data.rho_C_LR +*/
+ S_L * phase_transition_data.drho_C_LR_dp_GR;
+
+ dfC_2a.dp_GR = -porosity_data.phi * phase_transition_data.drho_C_GR_dp_GR -
+ drho_C_FR_dp_GR * pCap() *
+ (biot_data() - porosity_data.phi) * beta_p_SR();
+
+ double const dS_G_dp_cap = -dS_L_dp_cap();
+ double const rho_C_FR = S_G * constituent_density_data.rho_C_GR +
+ S_L * constituent_density_data.rho_C_LR;
+
+ // TODO (naumov) Extend for partially saturated media.
+ constexpr double drho_C_GR_dp_cap = 0;
+
+ double const drho_C_FR_dp_cap =
+ dS_G_dp_cap * constituent_density_data.rho_C_GR +
+ S_G * drho_C_GR_dp_cap +
+ dS_L_dp_cap() * constituent_density_data.rho_C_LR -
+ S_L * phase_transition_data.drho_C_LR_dp_LR;
+
+ dfC_2a.dp_cap =
+ porosity_data.phi *
+ (-phase_transition_data.drho_C_LR_dp_LR - drho_C_GR_dp_cap) -
+ drho_C_FR_dp_cap * pCap() * (biot_data() - porosity_data.phi) *
+ beta_p_SR() +
+ rho_C_FR * (biot_data() - porosity_data.phi) * beta_p_SR();
+
+ double const drho_C_FR_dT = S_G * phase_transition_data.drho_C_GR_dT +
+ S_L * phase_transition_data.drho_C_LR_dT;
+ dfC_2a.dT = porosity_d_data.dphi_dT * (constituent_density_data.rho_C_LR -
+ constituent_density_data.rho_C_GR) +
+ porosity_data.phi * (phase_transition_data.drho_C_LR_dT -
+ phase_transition_data.drho_C_GR_dT) -
+ drho_C_FR_dT * pCap() * (biot_data() - porosity_data.phi) *
+ beta_p_SR() +
+ rho_C_FR * pCap() * porosity_d_data.dphi_dT * beta_p_SR();
+}
+
+void FC3aModel::eval(
+ double const dt,
+ ConstituentDensityData const& constituent_density_data,
+ PrevState<ConstituentDensityData> const& constituent_density_data_prev,
+ SaturationData const& S_L_data,
+ FC3aData& fC_3a) const
+{
+ if (dt == 0.)
+ {
+ fC_3a.a = 0;
+ return;
+ }
+
+ double const rho_C_GR_dot = (constituent_density_data.rho_C_GR -
+ constituent_density_data_prev->rho_C_GR) /
+ dt;
+ double const rho_C_LR_dot = (constituent_density_data.rho_C_LR -
+ constituent_density_data_prev->rho_C_LR) /
+ dt;
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ fC_3a.a = S_G * rho_C_GR_dot + S_L * rho_C_LR_dot;
+}
+
+void FC3aModel::dEval(
+ double const dt,
+ ConstituentDensityData const& constituent_density_data,
+ PrevState<ConstituentDensityData> const& constituent_density_data_prev,
+ PhaseTransitionData const& phase_transition_data,
+ SaturationData const& S_L_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ FC3aDerivativeData& dfC_3a) const
+{
+ if (dt == 0.)
+ {
+ dfC_3a.dp_GR = 0.;
+ dfC_3a.dp_cap = 0.;
+ dfC_3a.dT = 0.;
+ return;
+ }
+ double const rho_C_GR_dot = (constituent_density_data.rho_C_GR -
+ constituent_density_data_prev->rho_C_GR) /
+ dt;
+ double const rho_C_LR_dot = (constituent_density_data.rho_C_LR -
+ constituent_density_data_prev->rho_C_LR) /
+ dt;
+
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ dfC_3a.dp_GR =
+ /*(dS_G_dp_GR = 0) * rho_C_GR_dot +*/
+ S_G * phase_transition_data.drho_C_GR_dp_GR / dt +
+ /*(dS_L_dp_GR = 0) * rho_C_LR_dot +*/
+ S_L * phase_transition_data.drho_C_LR_dp_GR / dt;
+
+ double const dS_G_dp_cap = -dS_L_dp_cap();
+ // TODO (naumov) Extend for partially saturated media.
+ constexpr double drho_C_GR_dp_cap = 0;
+
+ dfC_3a.dp_cap = dS_G_dp_cap * rho_C_GR_dot + S_G * drho_C_GR_dp_cap / dt +
+ dS_L_dp_cap() * rho_C_LR_dot -
+ S_L * phase_transition_data.drho_C_LR_dp_LR / dt;
+
+ dfC_3a.dT = S_G * phase_transition_data.drho_C_GR_dT / dt +
+ S_L * phase_transition_data.drho_C_LR_dT / dt;
+}
+
+template <int DisplacementDim>
+void FC4LCpGModel<DisplacementDim>::eval(
+ AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FC4LCpGData<DisplacementDim>& fC_4_LCpG) const
+{
+ GlobalDimMatrix<DisplacementDim> const advection_C =
+ advection_data.advection_C_G + advection_data.advection_C_L;
+
+ double const sD_G = phase_transition_data.diffusion_coefficient_vapour;
+ double const sD_L = phase_transition_data.diffusion_coefficient_solute;
+
+ double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi;
+ double const phi_L = S_L_data.S_L * porosity_data.phi;
+
+ double const diffusion_CGpGR = -phi_G * fluid_density_data.rho_GR * sD_G *
+ phase_transition_data.dxmWG_dpGR;
+ double const diffusion_CLpGR = -phi_L * fluid_density_data.rho_LR * sD_L *
+ phase_transition_data.dxmWL_dpGR;
+
+ double const diffusion_C_pGR = diffusion_CGpGR + diffusion_CLpGR;
+
+ auto const I =
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
+ fC_4_LCpG.L.noalias() = diffusion_C_pGR * I + advection_C;
+}
+
+template <int DisplacementDim>
+void FC4LCpGModel<DisplacementDim>::dEval(
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ ViscosityData const& viscosity_data,
+ PhaseTransitionData const& phase_transition_data,
+ AdvectionDerivativeData<DisplacementDim> const& advection_d_data,
+ FC4LCpGDerivativeData<DisplacementDim>& dfC_4_LCpG) const
+{
+ dfC_4_LCpG.dp_GR = advection_d_data.dadvection_C_dp_GR
+ // + ddiffusion_C_p_dp_GR TODO (naumov)
+ ;
+
+ dfC_4_LCpG.dp_cap = advection_d_data.dadvection_C_dp_cap
+ // + ddiffusion_C_p_dp_cap TODO (naumov)
+ ;
+
+ GlobalDimMatrix<DisplacementDim> const k_over_mu_G =
+ permeability_data.Ki * permeability_data.k_rel_G / viscosity_data.mu_GR;
+ GlobalDimMatrix<DisplacementDim> const k_over_mu_L =
+ permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR;
+
+ dfC_4_LCpG.dT = phase_transition_data.drho_C_GR_dT * k_over_mu_G +
+ phase_transition_data.drho_C_LR_dT * k_over_mu_L
+ // + ddiffusion_C_p_dT TODO (naumov)
+ ;
+}
+
+template struct FC4LCpGModel<2>;
+template struct FC4LCpGModel<3>;
+
+template <int DisplacementDim>
+void FC4LCpCModel<DisplacementDim>::eval(
+ AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FC4LCpCData<DisplacementDim>& fC_4_LCpC) const
+{
+ double const sD_G = phase_transition_data.diffusion_coefficient_vapour;
+ double const sD_L = phase_transition_data.diffusion_coefficient_solute;
+
+ double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi;
+ double const phi_L = S_L_data.S_L * porosity_data.phi;
+
+ double const diffusion_CGpCap = -phi_G * fluid_density_data.rho_GR * sD_G *
+ phase_transition_data.dxmWG_dpCap;
+ double const diffusion_CLpCap = -phi_L * fluid_density_data.rho_LR * sD_L *
+ phase_transition_data.dxmWL_dpCap;
+
+ double const diffusion_C_pCap = diffusion_CGpCap + diffusion_CLpCap;
+
+ auto const I =
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
+
+ fC_4_LCpC.L.noalias() = diffusion_C_pCap * I - advection_data.advection_C_L;
+}
+
+template <int DisplacementDim>
+void FC4LCpCModel<DisplacementDim>::dEval(
+ ConstituentDensityData const& constituent_density_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ PhaseTransitionData const& phase_transition_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ ViscosityData const& viscosity_data,
+ FC4LCpCDerivativeData<DisplacementDim>& dfC_4_LCpC) const
+{
+ ////// Diffusion Part /////
+ // TODO (naumov) d(diffusion_C_pCap)/dX for dxmW*/d* != 0
+
+ ////// Advection part /////
+ GlobalDimMatrix<DisplacementDim> const k_over_mu_L =
+ permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR;
+
+ dfC_4_LCpC.dp_GR = phase_transition_data.drho_C_LR_dp_GR * k_over_mu_L
+ //+ rhoCLR * (dk_over_mu_L_dp_GR = 0)
+ ;
+
+ auto const dk_over_mu_L_dp_cap = permeability_data.Ki *
+ permeability_data.dk_rel_L_dS_L *
+ dS_L_dp_cap() / viscosity_data.mu_LR;
+
+ dfC_4_LCpC.dp_cap = -phase_transition_data.drho_C_LR_dp_LR * k_over_mu_L +
+ constituent_density_data.rho_C_LR * dk_over_mu_L_dp_cap;
+
+ dfC_4_LCpC.dT = phase_transition_data.drho_W_LR_dT * k_over_mu_L
+ //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T))
+ ;
+}
+
+template struct FC4LCpCModel<2>;
+template struct FC4LCpCModel<3>;
+
+template <int DisplacementDim>
+void FC4LCTModel<DisplacementDim>::eval(
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FC4LCTData<DisplacementDim>& fC_4_LCT) const
+{
+ double const sD_G = phase_transition_data.diffusion_coefficient_vapour;
+ double const sD_L = phase_transition_data.diffusion_coefficient_solute;
+
+ double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi;
+ double const phi_L = S_L_data.S_L * porosity_data.phi;
+
+ double const diffusion_C_G_T = -phi_G * fluid_density_data.rho_GR * sD_G *
+ phase_transition_data.dxmWG_dT;
+ double const diffusion_C_L_T = -phi_L * fluid_density_data.rho_LR * sD_L *
+ phase_transition_data.dxmWL_dT;
+
+ double const diffusion_C_T = diffusion_C_G_T + diffusion_C_L_T;
+
+ auto const I =
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
+
+ fC_4_LCT.L.noalias() = diffusion_C_T * I;
+}
+
+template struct FC4LCTModel<2>;
+template struct FC4LCTModel<3>;
+
+void FC4MCpGModel::eval(BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const& beta_p_SR,
+ FC4MCpGData& fC_4_MCpG) const
+{
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ double const rho_C_FR = S_G * constituent_density_data.rho_C_GR +
+ S_L * constituent_density_data.rho_C_LR;
+
+ fC_4_MCpG.m = rho_C_FR * (biot_data() - porosity_data.phi) * beta_p_SR();
+}
+
+void FC4MCpGModel::dEval(BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PorosityDerivativeData const& porosity_d_data,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const& beta_p_SR,
+ FC4MCpGDerivativeData& dfC_4_MCpG) const
+{
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ double const rho_C_FR = S_G * constituent_density_data.rho_C_GR +
+ S_L * constituent_density_data.rho_C_LR;
+
+ double const drho_C_FR_dp_GR =
+ /*(dS_G_dp_GR = 0) * constituent_density_data.rho_C_GR +*/
+ S_G * phase_transition_data.drho_C_GR_dp_GR +
+ /*(dS_L_dp_GR = 0) * constituent_density_data.rho_C_LR +*/
+ S_L * phase_transition_data.drho_C_LR_dp_GR;
+
+ dfC_4_MCpG.dp_GR =
+ drho_C_FR_dp_GR * (biot_data() - porosity_data.phi) * beta_p_SR();
+
+ double const drho_C_FR_dT = S_G * phase_transition_data.drho_C_GR_dT +
+ S_L * phase_transition_data.drho_C_LR_dT;
+ dfC_4_MCpG.dT =
+ drho_C_FR_dT * (biot_data() - porosity_data.phi) * beta_p_SR() -
+ rho_C_FR * porosity_d_data.dphi_dT * beta_p_SR();
+}
+
+void FC4MCpCModel::eval(BiotData const& biot_data,
+ CapillaryPressureData const pCap,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ PrevState<SaturationData> const& S_L_data_prev,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const& beta_p_SR,
+ FC4MCpCData& fC_4_MCpC) const
+{
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ double const rho_C_FR = S_G * constituent_density_data.rho_C_GR +
+ S_L * constituent_density_data.rho_C_LR;
+
+ fC_4_MCpC.m =
+ -rho_C_FR * (biot_data() - porosity_data.phi) * beta_p_SR() * S_L;
+
+ fC_4_MCpC.ml =
+ (porosity_data.phi * (constituent_density_data.rho_C_LR -
+ constituent_density_data.rho_C_GR) -
+ rho_C_FR * pCap() * (biot_data() - porosity_data.phi) * beta_p_SR()) *
+ (S_L - S_L_data_prev->S_L);
+}
+
+template <int DisplacementDim>
+void FC4MCTModel<DisplacementDim>::eval(
+ BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ FC4MCTData& fC_4_MCT) const
+{
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ double const rho_C_FR = S_G * constituent_density_data.rho_C_GR +
+ S_L * constituent_density_data.rho_C_LR;
+
+ fC_4_MCT.m = -rho_C_FR * (biot_data() - porosity_data.phi) *
+ s_therm_exp_data.beta_T_SR;
+}
+
+template <int DisplacementDim>
+void FC4MCTModel<DisplacementDim>::dEval(
+ BiotData const& biot_data,
+ [[maybe_unused]] ConstituentDensityData const& constituent_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ [[maybe_unused]] PorosityDerivativeData const& porosity_d_data,
+ SaturationData const& S_L_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ FC4MCTDerivativeData& dfC_4_MCT) const
+{
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
+ double const rho_C_FR = S_G * constituent_density_data.rho_C_GR +
+ S_L * constituent_density_data.rho_C_LR;
+#endif
+
+ double const drho_C_FR_dT = S_G * phase_transition_data.drho_C_GR_dT +
+ S_L * phase_transition_data.drho_C_LR_dT;
+
+ dfC_4_MCT.dT = drho_C_FR_dT * (biot_data() - porosity_data.phi) *
+ s_therm_exp_data.beta_T_SR
+#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
+ + rho_C_FR * (biot_data() - porosity_d_data.dphi_dT) *
+ s_therm_exp_data.beta_T_SR
+#endif
+ ;
+}
+
+template struct FC4MCTModel<2>;
+template struct FC4MCTModel<3>;
+
+void FC4MCuModel::eval(BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ SaturationData const& S_L_data,
+ FC4MCuData& fC_4_MCu) const
+{
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ double const rho_C_FR = S_G * constituent_density_data.rho_C_GR +
+ S_L * constituent_density_data.rho_C_LR;
+
+ fC_4_MCu.m = rho_C_FR * biot_data();
+}
+
+void FC4MCuModel::dEval(BiotData const& biot_data,
+ PhaseTransitionData const& phase_transition_data,
+ SaturationData const& S_L_data,
+ FC4MCuDerivativeData& dfC_4_MCu) const
+{
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ double const drho_C_FR_dT = S_G * phase_transition_data.drho_C_GR_dT +
+ S_L * phase_transition_data.drho_C_LR_dT;
+ dfC_4_MCu.dT = drho_C_FR_dT * biot_data();
+}
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/CEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/CEquation.h
new file mode 100644
index 0000000000000000000000000000000000000000..c4de9c36a6402ba9338fcccfa908e7d8b7c64ef1
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/CEquation.h
@@ -0,0 +1,301 @@
+/**
+ * \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 "Advection.h"
+#include "Base.h"
+#include "Biot.h"
+#include "ConstitutiveDensity.h"
+#include "FluidDensity.h"
+#include "PermeabilityData.h"
+#include "PhaseTransitionData.h"
+#include "Porosity.h"
+#include "Saturation.h"
+#include "SolidCompressibility.h"
+#include "Viscosity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+struct FC1Data
+{
+ GlobalDimMatrix<DisplacementDim> A;
+};
+
+template <int DisplacementDim>
+struct FC1Model
+{
+ void eval(AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ FC1Data<DisplacementDim>& fC_1) const;
+};
+
+extern template struct FC1Model<2>;
+extern template struct FC1Model<3>;
+
+struct FC2aData
+{
+ double a = nan;
+};
+
+struct FC2aDerivativeData
+{
+ double dp_GR = nan;
+ double dp_cap = nan;
+ double dT = nan;
+};
+
+struct FC2aModel
+{
+ void eval(BiotData const biot_data,
+ CapillaryPressureData const pCap,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const beta_p_SR,
+ FC2aData& fC_2a) const;
+
+ void dEval(BiotData const& biot_data,
+ CapillaryPressureData const pCap,
+ ConstituentDensityData const& constituent_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PorosityDerivativeData const& porosity_d_data,
+ SaturationData const& S_L_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ SolidCompressibilityData const& beta_p_SR,
+ FC2aDerivativeData& dfC_2a) const;
+};
+
+struct FC3aData
+{
+ double a = nan;
+};
+
+struct FC3aDerivativeData
+{
+ double dp_GR = nan;
+ double dp_cap = nan;
+ double dT = nan;
+};
+
+struct FC3aModel
+{
+ void eval(
+ double const dt,
+ ConstituentDensityData const& constituent_density_data,
+ PrevState<ConstituentDensityData> const& constituent_density_data_prev,
+ SaturationData const& S_L_data,
+ FC3aData& fC_3a) const;
+
+ void dEval(
+ double const dt,
+ ConstituentDensityData const& constituent_density_data,
+ PrevState<ConstituentDensityData> const& constituent_density_data_prev,
+ PhaseTransitionData const& phase_transition_data,
+ SaturationData const& S_L_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ FC3aDerivativeData& dfC_3a) const;
+};
+
+template <int DisplacementDim>
+struct FC4LCpGData
+{
+ GlobalDimMatrix<DisplacementDim> L;
+};
+
+template <int DisplacementDim>
+struct FC4LCpGDerivativeData
+{
+ GlobalDimMatrix<DisplacementDim> dp_GR;
+ GlobalDimMatrix<DisplacementDim> dp_cap;
+ GlobalDimMatrix<DisplacementDim> dT;
+};
+
+template <int DisplacementDim>
+struct FC4LCpGModel
+{
+ void eval(AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FC4LCpGData<DisplacementDim>& fC_4_LCpG) const;
+
+ void dEval(PermeabilityData<DisplacementDim> const& permeability_data,
+ ViscosityData const& viscosity_data,
+ PhaseTransitionData const& phase_transition_data,
+ AdvectionDerivativeData<DisplacementDim> const& advection_d_data,
+ FC4LCpGDerivativeData<DisplacementDim>& dfC_4_LCpG) const;
+};
+
+extern template struct FC4LCpGModel<2>;
+extern template struct FC4LCpGModel<3>;
+
+template <int DisplacementDim>
+struct FC4LCpCData
+{
+ GlobalDimMatrix<DisplacementDim> L;
+};
+
+template <int DisplacementDim>
+struct FC4LCpCDerivativeData
+{
+ GlobalDimMatrix<DisplacementDim> dp_GR;
+ GlobalDimMatrix<DisplacementDim> dp_cap;
+ GlobalDimMatrix<DisplacementDim> dT;
+};
+
+template <int DisplacementDim>
+struct FC4LCpCModel
+{
+ void eval(AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FC4LCpCData<DisplacementDim>& fC_4_LCpC) const;
+
+ void dEval(ConstituentDensityData const& constituent_density_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ PhaseTransitionData const& phase_transition_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ ViscosityData const& viscosity_data,
+ FC4LCpCDerivativeData<DisplacementDim>& dfC_4_LCpC) const;
+};
+
+extern template struct FC4LCpCModel<2>;
+extern template struct FC4LCpCModel<3>;
+
+template <int DisplacementDim>
+struct FC4LCTData
+{
+ GlobalDimMatrix<DisplacementDim> L;
+};
+
+template <int DisplacementDim>
+struct FC4LCTModel
+{
+ void eval(FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FC4LCTData<DisplacementDim>& fC_4_LCT) const;
+};
+
+struct FC4MCpGData
+{
+ double m = nan;
+};
+
+struct FC4MCpGDerivativeData
+{
+ double dp_GR = nan;
+ double dT = nan;
+};
+
+struct FC4MCpGModel
+{
+ void eval(BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const& beta_p_SR,
+ FC4MCpGData& fC_4_MCpG) const;
+
+ void dEval(BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PorosityDerivativeData const& porosity_d_data,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const& beta_p_SR,
+ FC4MCpGDerivativeData& dfC_4_MCpG) const;
+};
+
+struct FC4MCpCData
+{
+ double m = nan;
+ double ml = nan;
+};
+
+struct FC4MCpCModel
+{
+ void eval(BiotData const& biot_data,
+ CapillaryPressureData const pCap,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ PrevState<SaturationData> const& S_L_data_prev,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const& beta_p_SR,
+ FC4MCpCData& fC_4_MCpC) const;
+};
+
+struct FC4MCTData
+{
+ double m = nan;
+};
+
+struct FC4MCTDerivativeData
+{
+ double dT = nan;
+};
+
+template <int DisplacementDim>
+struct FC4MCTModel
+{
+ void eval(
+ BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ FC4MCTData& fC_4_MCT) const;
+
+ void dEval(
+ BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PorosityDerivativeData const& porosity_d_data,
+ SaturationData const& S_L_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ FC4MCTDerivativeData& dfC_4_MCT) const;
+};
+
+struct FC4MCuData
+{
+ double m = nan;
+};
+
+struct FC4MCuDerivativeData
+{
+ double dT = nan;
+};
+
+struct FC4MCuModel
+{
+ void eval(BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ SaturationData const& S_L_data,
+ FC4MCuData& fC_4_MCu) const;
+
+ void dEval(BiotData const& biot_data,
+ PhaseTransitionData const& phase_transition_data,
+ SaturationData const& S_L_data,
+ FC4MCuDerivativeData& dfC_4_MCu) const;
+};
+
+extern template struct FC4MCTModel<2>;
+extern template struct FC4MCTModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index 7c4611cea705da0a547800e0e7999bc1fd826a7c..38f09a600852914be665fefb871965609833f55c 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -9,8 +9,10 @@
#pragma once
+#include "Advection.h"
#include "Biot.h"
#include "Bishops.h"
+#include "CEquation.h"
#include "ConstitutiveDensity.h"
#include "DarcyVelocity.h"
#include "DiffusionVelocity.h"
@@ -18,6 +20,7 @@
#include "Enthalpy.h"
#include "EquivalentPlasticStrainData.h"
#include "FluidDensity.h"
+#include "Gravity.h"
#include "InternalEnergy.h"
#include "MassMoleFractions.h"
#include "MechanicalStrain.h"
@@ -35,10 +38,13 @@
#include "SolidMechanics.h"
#include "SolidThermalExpansion.h"
#include "Swelling.h"
+#include "TEquation.h"
#include "ThermalConductivity.h"
#include "TotalStress.h"
+#include "UEquation.h"
#include "VapourPartialPressure.h"
#include "Viscosity.h"
+#include "WEquation.h"
namespace ProcessLib::TH2M
{
@@ -99,7 +105,8 @@ struct OutputData
{
ProcessLib::ConstitutiveRelations::StrainData<DisplacementDim> eps_data;
PermeabilityData<DisplacementDim> permeability_data;
- EnthalpyData enthalpy_data;
+ FluidEnthalpyData fluid_enthalpy_data;
+ SolidEnthalpyData solid_enthalpy_data;
MassMoleFractionsData mass_mole_fractions_data;
FluidDensityData fluid_density_data;
VapourPartialPressureData vapour_pressure_data;
@@ -114,7 +121,8 @@ struct OutputData
return Reflection::reflectWithoutName(&Self::eps_data,
&Self::permeability_data,
- &Self::enthalpy_data,
+ &Self::fluid_enthalpy_data,
+ &Self::solid_enthalpy_data,
&Self::mass_mole_fractions_data,
&Self::fluid_density_data,
&Self::vapour_pressure_data,
@@ -141,73 +149,76 @@ struct ConstitutiveTempData
ElasticTangentStiffnessData<DisplacementDim> C_el_data;
BiotData biot_data;
SolidCompressibilityData beta_p_SR;
- SaturationDataDeriv dS_L_dp_cap;
BishopsData chi_S_L;
SolidThermalExpansionData<DisplacementDim> s_therm_exp_data;
TotalStressData<DisplacementDim> total_stress_data;
EquivalentPlasticStrainData equivalent_plastic_strain_data;
ViscosityData viscosity_data;
PhaseTransitionData phase_transition_data;
- PorosityDerivativeData porosity_d_data;
- SolidDensityDerivativeData solid_density_d_data;
SolidHeatCapacityData solid_heat_capacity_data;
ThermalConductivityData<DisplacementDim> thermal_conductivity_data;
EffectiveVolumetricEnthalpy effective_volumetric_enthalpy_data;
- EffectiveVolumetricEnthalpyDerivatives effective_volumetric_enthalpy_d_data;
+ AdvectionData<DisplacementDim> advection_data;
+ VolumetricBodyForce<DisplacementDim> volumetric_body_force;
+ FC1Data<DisplacementDim> fC_1;
+ FC2aData fC_2a;
+ FC3aData fC_3a;
+ FC4LCpGData<DisplacementDim> fC_4_LCpG;
+ FC4LCpCData<DisplacementDim> fC_4_LCpC;
+ FC4LCTData<DisplacementDim> fC_4_LCT;
+ FC4MCpGData fC_4_MCpG;
+ FC4MCpCData fC_4_MCpC;
+ FC4MCTData fC_4_MCT;
+ FC4MCuData fC_4_MCu;
+
+ FW1Data<DisplacementDim> fW_1;
+ FW2Data fW_2;
+ FW3aData fW_3a;
+ FW4LWpGData<DisplacementDim> fW_4_LWpG;
+ FW4LWpCData<DisplacementDim> fW_4_LWpC;
+ FW4LWTData<DisplacementDim> fW_4_LWT;
+ FW4MWpGData fW_4_MWpG;
+ FW4MWpCData fW_4_MWpC;
+ FW4MWTData fW_4_MWT;
+ FW4MWuData fW_4_MWu;
+
+ FT1Data fT_1;
+ FT2Data<DisplacementDim> fT_2;
+ FT3Data<DisplacementDim> fT_3;
+
+ FU2KUpCData fu_2_KupC;
+};
+
+/// Data that stores intermediate values (derivatives), which are not needed
+/// outside the constitutive setting.
+template <int DisplacementDim>
+struct DerivativesData
+{
+ SaturationDataDeriv dS_L_dp_cap;
+ AdvectionDerivativeData<DisplacementDim> advection_d_data;
+ PorosityDerivativeData porosity_d_data;
+ ThermalConductivityDerivativeData<DisplacementDim>
+ thermal_conductivity_d_data;
+ SolidDensityDerivativeData solid_density_d_data;
EffectiveVolumetricInternalEnergyDerivatives
effective_volumetric_internal_energy_d_data;
-
- using DisplacementDimVector = Eigen::Matrix<double, DisplacementDim, 1>;
- using DisplacementDimMatrix =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
-
- DisplacementDimVector drho_GR_h_w_eff_dp_GR_Npart;
- DisplacementDimMatrix drho_GR_h_w_eff_dp_GR_gradNpart;
- DisplacementDimVector drho_LR_h_w_eff_dp_cap_Npart;
- DisplacementDimMatrix drho_LR_h_w_eff_dp_cap_gradNpart;
- DisplacementDimVector drho_GR_h_w_eff_dT;
- DisplacementDimMatrix dfW_4_LWpG_a_dp_GR;
- DisplacementDimMatrix dfW_4_LWpG_a_dp_cap;
- DisplacementDimMatrix dfW_4_LWpG_a_dT;
- DisplacementDimMatrix dfW_4_LWpG_d_dp_GR;
- DisplacementDimMatrix dfW_4_LWpG_d_dp_cap;
- DisplacementDimMatrix dfW_4_LWpG_d_dT;
- DisplacementDimMatrix dfW_4_LWpC_a_dp_GR;
- DisplacementDimMatrix dfW_4_LWpC_a_dp_cap;
- DisplacementDimMatrix dfW_4_LWpC_a_dT;
- DisplacementDimMatrix dfW_4_LWpC_d_dp_GR;
- DisplacementDimMatrix dfW_4_LWpC_d_dp_cap;
- DisplacementDimMatrix dfW_4_LWpC_d_dT;
- DisplacementDimMatrix dfC_4_LCpG_dT;
- DisplacementDimMatrix dfC_4_LCpC_a_dp_GR;
- DisplacementDimMatrix dfC_4_LCpC_a_dp_cap;
- DisplacementDimMatrix dfC_4_LCpC_a_dT;
- DisplacementDimMatrix dfC_4_LCpC_d_dp_GR;
- DisplacementDimMatrix dfC_4_LCpC_d_dp_cap;
- DisplacementDimMatrix dfC_4_LCpC_d_dT;
- DisplacementDimMatrix dadvection_C_dp_GR;
- DisplacementDimMatrix dadvection_C_dp_cap;
- DisplacementDimMatrix dk_over_mu_G_dp_cap;
- DisplacementDimMatrix dk_over_mu_L_dp_cap;
- double dfC_4_MCpG_dp_GR = std::numeric_limits<double>::quiet_NaN();
- double dfC_4_MCpG_dT = std::numeric_limits<double>::quiet_NaN();
- double dfC_4_MCT_dT = std::numeric_limits<double>::quiet_NaN();
- double dfC_4_MCu_dT = std::numeric_limits<double>::quiet_NaN();
- double dfC_3a_dp_GR = std::numeric_limits<double>::quiet_NaN();
- double dfC_3a_dp_cap = std::numeric_limits<double>::quiet_NaN();
- double dfC_3a_dT = std::numeric_limits<double>::quiet_NaN();
- double dfC_2a_dp_GR = std::numeric_limits<double>::quiet_NaN();
- double dfC_2a_dp_cap = std::numeric_limits<double>::quiet_NaN();
- double dfC_2a_dT = std::numeric_limits<double>::quiet_NaN();
- double dfW_2a_dp_GR = std::numeric_limits<double>::quiet_NaN();
- double dfW_2b_dp_GR = std::numeric_limits<double>::quiet_NaN();
- double dfW_2a_dp_cap = std::numeric_limits<double>::quiet_NaN();
- double dfW_2b_dp_cap = std::numeric_limits<double>::quiet_NaN();
- double dfW_2a_dT = std::numeric_limits<double>::quiet_NaN();
- double dfW_2b_dT = std::numeric_limits<double>::quiet_NaN();
- double dfW_3a_dp_GR = std::numeric_limits<double>::quiet_NaN();
- double dfW_3a_dp_cap = std::numeric_limits<double>::quiet_NaN();
- double dfW_3a_dT = std::numeric_limits<double>::quiet_NaN();
+ EffectiveVolumetricEnthalpyDerivatives effective_volumetric_enthalpy_d_data;
+ FC2aDerivativeData dfC_2a;
+ FC3aDerivativeData dfC_3a;
+ FC4LCpGDerivativeData<DisplacementDim> dfC_4_LCpG;
+ FC4LCpCDerivativeData<DisplacementDim> dfC_4_LCpC;
+ FC4MCpGDerivativeData dfC_4_MCpG;
+ FC4MCTDerivativeData dfC_4_MCT;
+ FC4MCuDerivativeData dfC_4_MCu;
+ FW2DerivativeData dfW_2;
+ FW3aDerivativeData dfW_3a;
+ FW4LWpGDerivativeData<DisplacementDim> dfW_4_LWpG;
+ FW4LWpCDerivativeData<DisplacementDim> dfW_4_LWpC;
+ FT1DerivativeData dfT_1;
+ FT2DerivativeData<DisplacementDim> dfT_2;
+ FU1KUTDerivativeData<DisplacementDim> dfu_1_KuT;
+ FU2KUpCDerivativeData dfu_2_KupC;
};
+
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index e37d7c906d204a71574007e1dea7db3a946e4bf4..277093d023f50498109c4f8c4955b90b9fb4cf79 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -9,9 +9,15 @@
#pragma once
+#include "Advection.h"
#include "Biot.h"
#include "Bishops.h"
+#include "CEquation.h"
+#include "DiffusionVelocity.h"
#include "ElasticTangentStiffnessModel.h"
+#include "Enthalpy.h"
+#include "Gravity.h"
+#include "InternalEnergy.h"
#include "MechanicalStrain.h"
#include "PermeabilityModel.h"
#include "PhaseTransitionModel.h"
@@ -24,9 +30,12 @@
#include "SolidMechanics.h"
#include "SolidThermalExpansion.h"
#include "Swelling.h"
+#include "TEquation.h"
#include "ThermalConductivity.h"
#include "TotalStress.h"
+#include "UEquation.h"
#include "Viscosity.h"
+#include "WEquation.h"
namespace ProcessLib::TH2M
{
@@ -74,6 +83,41 @@ struct ConstitutiveModels
#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
SolidHeatCapacityModel solid_heat_capacity_model;
ThermalConductivityModel<DisplacementDim> thermal_conductivity_model;
+ AdvectionModel<DisplacementDim> advection_model;
+ InternalEnergyModel internal_energy_model;
+ EffectiveVolumetricEnthalpyModel effective_volumetric_enthalpy_model;
+ GravityModel<DisplacementDim> gravity_model;
+ DiffusionVelocityModel<DisplacementDim> diffusion_velocity_model;
+ SolidEnthalpyModel solid_enthalpy_model;
+ DarcyVelocityModel<DisplacementDim> darcy_velocity_model;
+ FC1Model<DisplacementDim> fC_1_model;
+ FC2aModel fC_2a_model;
+ FC3aModel fC_3a_model;
+ FC4LCpGModel<DisplacementDim> fC_4_LCpG_model;
+ FC4LCpCModel<DisplacementDim> fC_4_LCpC_model;
+ FC4LCTModel<DisplacementDim> fC_4_LCT_model;
+ FC4MCpGModel fC_4_MCpG_model;
+ FC4MCpCModel fC_4_MCpC_model;
+ FC4MCTModel<DisplacementDim> fC_4_MCT_model;
+ FC4MCuModel fC_4_MCu_model;
+
+ FW1Model<DisplacementDim> fW_1_model;
+ FW2Model fW_2_model;
+ FW3aModel fW_3a_model;
+ FW4LWpGModel<DisplacementDim> fW_4_LWpG_model;
+ FW4LWpCModel<DisplacementDim> fW_4_LWpC_model;
+ FW4LWTModel<DisplacementDim> fW_4_LWT_model;
+ FW4MWpGModel fW_4_MWpG_model;
+ FW4MWpCModel fW_4_MWpC_model;
+ FW4MWTModel<DisplacementDim> fW_4_MWT_model;
+ FW4MWuModel fW_4_MWu_model;
+
+ FT1Model fT_1_model;
+ FT2Model<DisplacementDim> fT_2_model;
+ FT3Model<DisplacementDim> fT_3_model;
+
+ FU1KUTModel<DisplacementDim> fu_1_KuT_model;
+ FU2KUpCModel fu_2_KupC_model;
};
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5083d8c65c25e801bdcddc94ba51c3c2da8dcb59
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.cpp
@@ -0,0 +1,43 @@
+/**
+ * \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 "DarcyVelocity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void DarcyVelocityModel<DisplacementDim>::eval(
+ CapillaryPressureGradientData<DisplacementDim> const& grad_p_cap,
+ FluidDensityData const& fluid_density_data,
+ GasPressureGradientData<DisplacementDim> const& grad_p_GR,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+ ViscosityData const& viscosity_data,
+ DarcyVelocityData<DisplacementDim>& darcy_velocity_data) const
+{
+ auto const k_over_mu_G =
+ permeability_data.Ki * permeability_data.k_rel_G / viscosity_data.mu_GR;
+ auto const k_over_mu_L =
+ permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR;
+
+ darcy_velocity_data.w_GS =
+ k_over_mu_G * fluid_density_data.rho_GR * specific_body_force() -
+ k_over_mu_G * grad_p_GR();
+ darcy_velocity_data.w_LS =
+ k_over_mu_L * grad_p_cap() +
+ k_over_mu_L * fluid_density_data.rho_LR * specific_body_force() -
+ k_over_mu_L * grad_p_GR();
+}
+
+template struct DarcyVelocityModel<2>;
+template struct DarcyVelocityModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
index 84e66d31d02bcf744f6d5affa6a7183695d2392b..91751b830f485a309fea551dd463be2e9d7769d9 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
@@ -10,7 +10,10 @@
#pragma once
#include "Base.h"
+#include "FluidDensity.h"
+#include "PermeabilityData.h"
#include "ProcessLib/Reflection/ReflectionData.h"
+#include "Viscosity.h"
namespace ProcessLib::TH2M
{
@@ -32,5 +35,20 @@ struct DarcyVelocityData
R::makeReflectionData("velocity_liquid", &Self::w_LS)};
}
};
+
+template <int DisplacementDim>
+struct DarcyVelocityModel
+{
+ void eval(CapillaryPressureGradientData<DisplacementDim> const& grad_p_cap,
+ FluidDensityData const& fluid_density_data,
+ GasPressureGradientData<DisplacementDim> const& grad_p_GR,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+ ViscosityData const& viscosity_data,
+ DarcyVelocityData<DisplacementDim>& darcy_velocity_data) const;
+};
+
+extern template struct DarcyVelocityModel<2>;
+extern template struct DarcyVelocityModel<3>;
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..20bd591702698b46790543ed7d87a73b4512a6d0
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.cpp
@@ -0,0 +1,90 @@
+/**
+ * \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 "DiffusionVelocity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void DiffusionVelocityModel<DisplacementDim>::eval(
+ CapillaryPressureGradientData<DisplacementDim> const& grad_p_cap,
+ GasPressureGradientData<DisplacementDim> const& grad_p_GR,
+ MassMoleFractionsData const& mass_mole_fractions_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ TemperatureGradientData<DisplacementDim> const& grad_T,
+ DiffusionVelocityData<DisplacementDim>& diffusion_velocity_data) const
+{
+ auto const gradxmWG = phase_transition_data.dxmWG_dpGR * grad_p_GR() +
+ phase_transition_data.dxmWG_dpCap * grad_p_cap() +
+ phase_transition_data.dxmWG_dT * grad_T();
+ auto const gradxmCG = -gradxmWG;
+
+ auto const gradxmWL = phase_transition_data.dxmWL_dpGR * grad_p_GR() +
+ phase_transition_data.dxmWL_dpCap * grad_p_cap() +
+ phase_transition_data.dxmWL_dT * grad_T();
+ auto const gradxmCL = -gradxmWL;
+
+ double const phi_L = S_L_data.S_L * porosity_data.phi;
+ double const phi_G = (1. - S_L_data.S_L) * porosity_data.phi;
+
+ double const phi_G_D_vapour =
+ phi_G * phase_transition_data.diffusion_coefficient_vapour;
+ double const phi_L_D_solute =
+ phi_L * phase_transition_data.diffusion_coefficient_solute;
+
+ if (mass_mole_fractions_data.xmCG == 0)
+ {
+ diffusion_velocity_data.d_CG = GlobalDimVector<DisplacementDim>::Zero();
+ }
+ else
+ {
+ diffusion_velocity_data.d_CG =
+ -phi_G_D_vapour / mass_mole_fractions_data.xmCG * gradxmCG;
+ }
+
+ if (mass_mole_fractions_data.xmCG == 1)
+ {
+ diffusion_velocity_data.d_WG = GlobalDimVector<DisplacementDim>::Zero();
+ }
+ else
+ {
+ diffusion_velocity_data.d_WG =
+ -phi_G_D_vapour / (1 - mass_mole_fractions_data.xmCG) * gradxmWG;
+ }
+
+ if (mass_mole_fractions_data.xmWL == 1)
+
+ {
+ diffusion_velocity_data.d_CL = GlobalDimVector<DisplacementDim>::Zero();
+ }
+ else
+ {
+ diffusion_velocity_data.d_CL =
+ -phi_L_D_solute / (1. - mass_mole_fractions_data.xmWL) * gradxmCL;
+ }
+
+ if (mass_mole_fractions_data.xmWL == 0)
+ {
+ diffusion_velocity_data.d_WL = GlobalDimVector<DisplacementDim>::Zero();
+ }
+ else
+ {
+ diffusion_velocity_data.d_WL =
+ -phi_L_D_solute / mass_mole_fractions_data.xmWL * gradxmWL;
+ }
+}
+
+template struct DiffusionVelocityModel<2>;
+template struct DiffusionVelocityModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h
index 25a203ca2c0c03fee94a64527880c3e07f37929e..7ec45405ed82137de5b2d6ae8206c1ad10ff5833 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h
@@ -10,7 +10,11 @@
#pragma once
#include "Base.h"
+#include "MassMoleFractions.h"
+#include "PhaseTransitionData.h"
+#include "Porosity.h"
#include "ProcessLib/Reflection/ReflectionData.h"
+#include "Saturation.h"
namespace ProcessLib::TH2M
{
@@ -38,5 +42,22 @@ struct DiffusionVelocityData
&Self::d_WL)};
}
};
+
+template <int DisplacementDim>
+struct DiffusionVelocityModel
+{
+ void eval(
+ CapillaryPressureGradientData<DisplacementDim> const& grad_p_cap,
+ GasPressureGradientData<DisplacementDim> const& grad_p_GR,
+ MassMoleFractionsData const& mass_mole_fractions_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ TemperatureGradientData<DisplacementDim> const& grad_T,
+ DiffusionVelocityData<DisplacementDim>& diffusion_velocity_data) const;
+};
+
+extern template struct DiffusionVelocityModel<2>;
+extern template struct DiffusionVelocityModel<3>;
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..046bb5eb2fb37cadf600c599be808eae2b7ad4a4
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.cpp
@@ -0,0 +1,101 @@
+/**
+ * \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 "InternalEnergy.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+void EffectiveVolumetricEnthalpyModel::eval(
+ FluidDensityData const& fluid_density_data,
+ FluidEnthalpyData const& fluid_enthalpy_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidDensityData const& solid_density_data,
+ SolidEnthalpyData const& solid_enthalpy_data,
+ EffectiveVolumetricEnthalpy& effective_volumetric_enthalpy_data) const
+{
+ auto const phi_L = S_L_data.S_L * porosity_data.phi;
+ auto const phi_G = (1. - S_L_data.S_L) * porosity_data.phi;
+ double const phi_S = 1. - porosity_data.phi;
+
+ effective_volumetric_enthalpy_data.rho_h_eff =
+ phi_G * fluid_density_data.rho_GR * fluid_enthalpy_data.h_G +
+ phi_L * fluid_density_data.rho_LR * fluid_enthalpy_data.h_L +
+ phi_S * solid_density_data.rho_SR * solid_enthalpy_data.h_S;
+}
+
+void EffectiveVolumetricEnthalpyModel::dEval(
+ FluidDensityData const& fluid_density_data,
+ FluidEnthalpyData const& fluid_enthalpy_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PorosityDerivativeData const& porosity_d_data,
+ SaturationData const& S_L_data,
+ SolidDensityData const& solid_density_data,
+ SolidDensityDerivativeData const& solid_density_d_data,
+ SolidEnthalpyData const& solid_enthalpy_data,
+ SolidHeatCapacityData const& solid_heat_capacity_data,
+ EffectiveVolumetricEnthalpyDerivatives&
+ effective_volumetric_enthalpy_d_data) const
+{
+ auto const phi_L = S_L_data.S_L * porosity_data.phi;
+ auto const phi_G = (1. - S_L_data.S_L) * porosity_data.phi;
+ double const phi_S = 1. - porosity_data.phi;
+
+ // From p_LR = p_GR - p_cap it follows for
+ // drho_LR/dp_GR = drho_LR/dp_LR * dp_LR/dp_GR
+ // = drho_LR/dp_LR * (dp_GR/dp_GR - dp_cap/dp_GR)
+ // = drho_LR/dp_LR * (1 - 0)
+ double const drho_LR_dp_GR = phase_transition_data.drho_LR_dp_LR;
+ double const drho_LR_dp_cap = -phase_transition_data.drho_LR_dp_LR;
+ // drho_GR_dp_cap = 0;
+
+ effective_volumetric_enthalpy_d_data.drho_h_eff_dp_GR =
+ /*(dphi_G_dp_GR = 0) * fluid_density_data.rho_GR *
+ fluid_enthalpy_data.h_G +*/
+ phi_G * phase_transition_data.drho_GR_dp_GR * fluid_enthalpy_data.h_G +
+ /*(dphi_L_dp_GR = 0) * fluid_density_data.rho_LR *
+ fluid_enthalpy_data.h_L +*/
+ phi_L * drho_LR_dp_GR * fluid_enthalpy_data.h_L;
+ effective_volumetric_enthalpy_d_data.drho_h_eff_dp_cap =
+ porosity_d_data.dphi_L_dp_cap * fluid_density_data.rho_GR *
+ fluid_enthalpy_data.h_G +
+ /*phi_G * (drho_GR_dp_cap = 0) * fluid_enthalpy_data.h_G +*/
+ porosity_d_data.dphi_L_dp_cap * fluid_density_data.rho_LR *
+ fluid_enthalpy_data.h_L +
+ phi_L * drho_LR_dp_cap * fluid_enthalpy_data.h_L;
+
+ // TODO (naumov) Extend for temperature dependent porosities.
+ constexpr double dphi_G_dT = 0;
+ constexpr double dphi_L_dT = 0;
+ effective_volumetric_enthalpy_d_data.drho_h_eff_dT =
+ dphi_G_dT * fluid_density_data.rho_GR * fluid_enthalpy_data.h_G +
+ phi_G * phase_transition_data.drho_GR_dT * fluid_enthalpy_data.h_G +
+ phi_G * fluid_density_data.rho_GR * phase_transition_data.dh_G_dT +
+ dphi_L_dT * fluid_density_data.rho_LR * fluid_enthalpy_data.h_L +
+ phi_L * phase_transition_data.drho_LR_dT * fluid_enthalpy_data.h_L +
+ phi_L * fluid_density_data.rho_LR * phase_transition_data.dh_L_dT -
+ porosity_d_data.dphi_dT * solid_density_data.rho_SR *
+ solid_enthalpy_data.h_S +
+ phi_S * solid_density_d_data.drho_SR_dT * solid_enthalpy_data.h_S +
+ phi_S * solid_density_data.rho_SR * solid_heat_capacity_data();
+}
+
+void SolidEnthalpyModel::eval(
+ SolidHeatCapacityData const& solid_heat_capacity_data,
+ TemperatureData const& T_data,
+ SolidEnthalpyData& solid_enthalpy_data) const
+{
+ solid_enthalpy_data.h_S = solid_heat_capacity_data() * T_data.T;
+}
+
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h
index b7a1686618b93faa2386ab2df657d3ae91a92512..4bff673ce9abeb66f8566e4bc089f066a584a2cb 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h
@@ -10,26 +10,44 @@
#pragma once
#include "Base.h"
+#include "Enthalpy.h"
+#include "FluidDensity.h"
+#include "PhaseTransitionData.h"
+#include "Porosity.h"
#include "ProcessLib/Reflection/ReflectionData.h"
+#include "Saturation.h"
+#include "SolidDensity.h"
+#include "SolidHeatCapacity.h"
namespace ProcessLib::TH2M
{
namespace ConstitutiveRelations
{
-struct EnthalpyData
+struct FluidEnthalpyData
{
double h_G = nan;
double h_L = nan;
- double h_S = nan;
static auto reflect()
{
- using Self = EnthalpyData;
+ using Self = FluidEnthalpyData;
namespace R = ProcessLib::Reflection;
return std::tuple{R::makeReflectionData("enthalpy_gas", &Self::h_G),
- R::makeReflectionData("enthalpy_liquid", &Self::h_L),
- R::makeReflectionData("enthalpy_solid", &Self::h_S)};
+ R::makeReflectionData("enthalpy_liquid", &Self::h_L)};
+ }
+};
+
+struct SolidEnthalpyData
+{
+ double h_S = nan;
+
+ static auto reflect()
+ {
+ using Self = SolidEnthalpyData;
+ namespace R = ProcessLib::Reflection;
+
+ return std::tuple{R::makeReflectionData("enthalpy_solid", &Self::h_S)};
}
};
@@ -45,5 +63,36 @@ struct EffectiveVolumetricEnthalpyDerivatives
double drho_h_eff_dp_cap = nan;
};
+struct EffectiveVolumetricEnthalpyModel
+{
+ void eval(
+ FluidDensityData const& fluid_density_data,
+ FluidEnthalpyData const& fluid_enthalpy_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidDensityData const& solid_density_data,
+ SolidEnthalpyData const& solid_enthalpy_data,
+ EffectiveVolumetricEnthalpy& effective_volumetric_enthalpy_data) const;
+
+ void dEval(FluidDensityData const& fluid_density_data,
+ FluidEnthalpyData const& fluid_enthalpy_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PorosityDerivativeData const& porosity_d_data,
+ SaturationData const& S_L_data,
+ SolidDensityData const& solid_density_data,
+ SolidDensityDerivativeData const& solid_density_d_data,
+ SolidEnthalpyData const& solid_enthalpy_data,
+ SolidHeatCapacityData const& solid_heat_capacity_data,
+ EffectiveVolumetricEnthalpyDerivatives&
+ effective_volumetric_enthalpy_d_data) const;
+};
+
+struct SolidEnthalpyModel
+{
+ void eval(SolidHeatCapacityData const& solid_heat_capacity_data,
+ TemperatureData const& T_data,
+ SolidEnthalpyData& solid_enthalpy_data) const;
+};
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Gravity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Gravity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8a4e994a876afbe886e46993b78695c6a6ad455a
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Gravity.cpp
@@ -0,0 +1,39 @@
+/**
+ * \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 "Gravity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void GravityModel<DisplacementDim>::eval(
+ FluidDensityData const& fluid_density_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidDensityData const& solid_density_data,
+ SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+ VolumetricBodyForce<DisplacementDim>& volumetric_body_force) const
+{
+ auto const phi_G = (1. - S_L_data.S_L) * porosity_data.phi;
+ auto const phi_L = S_L_data.S_L * porosity_data.phi;
+ auto const phi_S = 1. - porosity_data.phi;
+
+ auto const rhoGR = fluid_density_data.rho_GR;
+ auto const rhoLR = fluid_density_data.rho_LR;
+ auto const rho =
+ phi_G * rhoGR + phi_L * rhoLR + phi_S * solid_density_data.rho_SR;
+ *volumetric_body_force = rho * specific_body_force();
+}
+
+template struct GravityModel<2>;
+template struct GravityModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Gravity.h b/ProcessLib/TH2M/ConstitutiveRelations/Gravity.h
new file mode 100644
index 0000000000000000000000000000000000000000..dc1aac245e673d7cb7d08f6a55bdde987b52a816
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Gravity.h
@@ -0,0 +1,42 @@
+/**
+ * \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"
+#include "FluidDensity.h"
+#include "Porosity.h"
+#include "Saturation.h"
+#include "SolidDensity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+using VolumetricBodyForce =
+ BaseLib::StrongType<GlobalDimVector<DisplacementDim>,
+ struct VolumetricBodyForceTag>;
+
+template <int DisplacementDim>
+struct GravityModel
+{
+ void eval(
+ FluidDensityData const& fluid_density_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidDensityData const& solid_density_data,
+ SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+ VolumetricBodyForce<DisplacementDim>& volumetric_body_force) const;
+};
+
+extern template struct GravityModel<2>;
+extern template struct GravityModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.cpp b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2aa2b1d75526d4454a24d8d54939b946b9505b46
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.cpp
@@ -0,0 +1,91 @@
+/**
+ * \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 "InternalEnergy.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+void InternalEnergyModel::eval(FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidDensityData const& solid_density_data,
+ SolidEnthalpyData const& solid_enthalpy_data,
+ InternalEnergyData& internal_energy_data) const
+{
+ auto const phi_L = S_L_data.S_L * porosity_data.phi;
+ auto const phi_G = (1. - S_L_data.S_L) * porosity_data.phi;
+ double const phi_S = 1. - porosity_data.phi;
+
+ auto const u_S = solid_enthalpy_data.h_S;
+
+ internal_energy_data() =
+ phi_G * fluid_density_data.rho_GR * phase_transition_data.uG +
+ phi_L * fluid_density_data.rho_LR * phase_transition_data.uL +
+ phi_S * solid_density_data.rho_SR * u_S;
+}
+
+void InternalEnergyModel::dEval(
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PorosityDerivativeData const& porosity_d_data,
+ SaturationData const& S_L_data,
+ SolidDensityData const& solid_density_data,
+ SolidDensityDerivativeData const& solid_density_d_data,
+ SolidEnthalpyData const& solid_enthalpy_data,
+ SolidHeatCapacityData const& solid_heat_capacity_data,
+ EffectiveVolumetricInternalEnergyDerivatives&
+ effective_volumetric_internal_energy_d_data) const
+{
+ auto const phi_L = S_L_data.S_L * porosity_data.phi;
+ auto const phi_G = (1. - S_L_data.S_L) * porosity_data.phi;
+ double const phi_S = 1. - porosity_data.phi;
+
+ auto const u_S = solid_enthalpy_data.h_S;
+ effective_volumetric_internal_energy_d_data.drho_u_eff_dT =
+ phi_G * phase_transition_data.drho_GR_dT * phase_transition_data.uG +
+ phi_G * fluid_density_data.rho_GR * phase_transition_data.du_G_dT +
+ phi_L * phase_transition_data.drho_LR_dT * phase_transition_data.uL +
+ phi_L * fluid_density_data.rho_LR * phase_transition_data.du_L_dT +
+ phi_S * solid_density_d_data.drho_SR_dT * u_S +
+ phi_S * solid_density_data.rho_SR * solid_heat_capacity_data() -
+ porosity_d_data.dphi_dT * solid_density_data.rho_SR * u_S;
+
+ // From p_LR = p_GR - p_cap it follows for
+ // drho_LR/dp_GR = drho_LR/dp_LR * dp_LR/dp_GR
+ // = drho_LR/dp_LR * (dp_GR/dp_GR - dp_cap/dp_GR)
+ // = drho_LR/dp_LR * (1 - 0)
+ double const drho_LR_dp_GR = phase_transition_data.drho_LR_dp_LR;
+ double const drho_LR_dp_cap = -phase_transition_data.drho_LR_dp_LR;
+ // drho_GR_dp_cap = 0;
+
+ effective_volumetric_internal_energy_d_data.drho_u_eff_dp_GR =
+ /*(dphi_G_dp_GR = 0) * fluid_density_data.rho_GR *
+ phase_transition_data.uG +*/
+ phi_G * phase_transition_data.drho_GR_dp_GR * phase_transition_data.uG +
+ phi_G * fluid_density_data.rho_GR * phase_transition_data.du_G_dp_GR +
+ /*(dphi_L_dp_GR = 0) * fluid_density_data.rho_LR *
+ phase_transition_data.uL +*/
+ phi_L * drho_LR_dp_GR * phase_transition_data.uL +
+ phi_L * fluid_density_data.rho_LR * phase_transition_data.du_L_dp_GR;
+
+ effective_volumetric_internal_energy_d_data.drho_u_eff_dp_cap =
+ -porosity_d_data.dphi_L_dp_cap * fluid_density_data.rho_GR *
+ phase_transition_data.uG +
+ /*phi_G * (drho_GR_dp_cap = 0) * phase_transition_data.uG +*/
+ porosity_d_data.dphi_L_dp_cap * fluid_density_data.rho_LR *
+ phase_transition_data.uL +
+ phi_L * drho_LR_dp_cap * phase_transition_data.uL +
+ phi_L * fluid_density_data.rho_LR * phase_transition_data.du_L_dp_cap;
+}
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h
index 73644146b3452dd53f22f0d4a3b3485d997b0e9a..d74dfd7167708c10cbe4a9beb146b30e7cabf3ec 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h
@@ -11,6 +11,13 @@
#include "Base.h"
#include "BaseLib/StrongType.h"
+#include "Enthalpy.h"
+#include "FluidDensity.h"
+#include "PhaseTransitionData.h"
+#include "Porosity.h"
+#include "Saturation.h"
+#include "SolidDensity.h"
+#include "SolidHeatCapacity.h"
namespace ProcessLib::TH2M
{
@@ -25,5 +32,28 @@ struct EffectiveVolumetricInternalEnergyDerivatives
double drho_u_eff_dp_GR = nan;
double drho_u_eff_dp_cap = nan;
};
+
+struct InternalEnergyModel
+{
+ void eval(FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidDensityData const& solid_density_data,
+ SolidEnthalpyData const& solid_enthalpy_data,
+ InternalEnergyData& internal_energy_data) const;
+
+ void dEval(FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PorosityDerivativeData const& porosity_d_data,
+ SaturationData const& S_L_data,
+ SolidDensityData const& solid_density_data,
+ SolidDensityDerivativeData const& solid_density_d_data,
+ SolidEnthalpyData const& solid_enthalpy_data,
+ SolidHeatCapacityData const& solid_heat_capacity_data,
+ EffectiveVolumetricInternalEnergyDerivatives&
+ effective_volumetric_internal_energy_d_data) const;
+};
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp
index 782b8b60ecbec18059c1bd9aef3134bd22c95baa..72e475c043427e0c843461edf2ef0e8ebe1987f9 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp
@@ -43,7 +43,7 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t,
CapillaryPressureData const& p_cap,
TemperatureData const& T_data,
PureLiquidDensityData const& rho_W_LR,
- EnthalpyData& enthalpy_data,
+ FluidEnthalpyData& fluid_enthalpy_data,
MassMoleFractionsData& mass_mole_fractions_data,
FluidDensityData& fluid_density_data,
VapourPartialPressureData& vapour_pressure_data,
@@ -108,14 +108,14 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t,
.template value<double>(variables, x_t.x, x_t.t, x_t.dt);
// specific phase enthalpies
- enthalpy_data.h_G = cpG * T;
- enthalpy_data.h_L = cpL * T;
+ fluid_enthalpy_data.h_G = cpG * T;
+ fluid_enthalpy_data.h_L = cpL * T;
cv.dh_G_dT = cpG;
cv.dh_L_dT = cpL;
// specific inner energies
- cv.uG = enthalpy_data.h_G - pGR / fluid_density_data.rho_GR;
- cv.uL = enthalpy_data.h_L;
+ cv.uG = fluid_enthalpy_data.h_G - pGR / fluid_density_data.rho_GR;
+ cv.uL = fluid_enthalpy_data.h_L;
cv.drho_GR_dT =
gas_phase[MaterialPropertyLib::PropertyType::density]
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h
index 918723e906012fb59a5bb1c1cd7c99b67dd23e0c..102443321b98a3df7afd77889eb6c04851c0ce4d 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h
@@ -28,7 +28,7 @@ struct NoPhaseTransition : PhaseTransitionModel
GasPressureData const& p_GR, CapillaryPressureData const& p_cap,
TemperatureData const& T_data,
PureLiquidDensityData const& rho_W_LR,
- EnthalpyData& enthalpy_data,
+ FluidEnthalpyData& fluid_enthalpy_data,
MassMoleFractionsData& mass_mole_fractions_data,
FluidDensityData& fluid_density_data,
VapourPartialPressureData& vapour_pressure_data,
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp
index 7f71fcb24920fee9f6fd6094cec574158f19a83d..597bdca0440849669630915b32c902fba17b7da0 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp
@@ -139,7 +139,7 @@ void PhaseTransition::eval(SpaceTimeData const& x_t,
CapillaryPressureData const& p_cap,
TemperatureData const& T_data,
PureLiquidDensityData const& rho_W_LR,
- EnthalpyData& enthalpy_data,
+ FluidEnthalpyData& fluid_enthalpy_data,
MassMoleFractionsData& mass_mole_fractions_data,
FluidDensityData& fluid_density_data,
VapourPartialPressureData& vapour_pressure_data,
@@ -333,11 +333,12 @@ void PhaseTransition::eval(SpaceTimeData const& x_t,
cv.hWG = cpWG * T + dh_evap;
// specific enthalpy of gas phase
- enthalpy_data.h_G = mass_mole_fractions_data.xmCG * cv.hCG + xmWG * cv.hWG;
+ fluid_enthalpy_data.h_G =
+ mass_mole_fractions_data.xmCG * cv.hCG + xmWG * cv.hWG;
cv.dh_G_dT = 0;
// specific inner energies of gas phase
- cv.uG = enthalpy_data.h_G - pGR / fluid_density_data.rho_GR;
+ cv.uG = fluid_enthalpy_data.h_G - pGR / fluid_density_data.rho_GR;
cv.du_G_dT = 0;
cv.du_G_dp_GR = 0;
@@ -505,11 +506,11 @@ void PhaseTransition::eval(SpaceTimeData const& x_t,
// specific enthalpy of liquid phase and its components
double const hCL = cpCL * T + dh_sol;
double const hWL = cpWL * T;
- enthalpy_data.h_L = xmCL * hCL + mass_mole_fractions_data.xmWL * hWL;
+ fluid_enthalpy_data.h_L = xmCL * hCL + mass_mole_fractions_data.xmWL * hWL;
cv.dh_L_dT = 0;
// specific inner energies of liquid phase
- cv.uL = enthalpy_data.h_L;
+ cv.uL = fluid_enthalpy_data.h_L;
cv.du_L_dT = 0;
// diffusion
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h
index 69d621b3dad66159bbf4b350d62e1102a74f893a..58766650f5347f837e947650a40ff257f435e069 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h
@@ -28,7 +28,7 @@ struct PhaseTransition : PhaseTransitionModel
GasPressureData const& p_GR, CapillaryPressureData const& p_cap,
TemperatureData const& T_data,
PureLiquidDensityData const& rho_W_LR,
- EnthalpyData& enthalpy_data,
+ FluidEnthalpyData& fluid_enthalpy_data,
MassMoleFractionsData& mass_mole_fractions_data,
FluidDensityData& fluid_density_data,
VapourPartialPressureData& vapour_pressure_data,
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h
index a0042c6afbbad7c47880838830b1ef95681387a5..b6e27f07574f0fc3b09227b5862add82d2748516 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h
@@ -51,10 +51,6 @@ struct PhaseTransitionData
double dxmWL_dpGR = nan;
double dxmWL_dpCap = nan;
double dxmWL_dT = nan;
- // TODO (naumov) This is zero in both models but used in the assembly.
- // Remove it and simplify assembly or correct the expressions in the phase
- // transition models?
- static constexpr double dxmWL_dpLR = 0.;
double diffusion_coefficient_vapour = nan;
double diffusion_coefficient_solute = nan;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h
index 8cf92e94f7e03421a82a9c6dd8f7c30178cfe37d..0642ce42ca85055a1142bca949a6191ae490d840 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h
@@ -52,7 +52,7 @@ struct PhaseTransitionModel
CapillaryPressureData const& p_cap,
TemperatureData const& T_data,
PureLiquidDensityData const& rho_W_LR,
- EnthalpyData& enthalpy_data,
+ FluidEnthalpyData& fluid_enthalpy_data,
MassMoleFractionsData& mass_mole_fractions_data,
FluidDensityData& fluid_density_data,
VapourPartialPressureData& vapour_pressure_data,
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp
index 508dc689d3a8c32a338a0cb77093b0a259c75a66..2688dd02698f3e4072a3a465d36a12badb4ee7e7 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp
@@ -16,8 +16,7 @@ namespace ConstitutiveRelations
void PorosityModel::eval(SpaceTimeData const& x_t,
MediaData const& media_data,
- PorosityData& porosity_data,
- PorosityDerivativeData& porosity_d_data) const
+ PorosityData& porosity_data) const
{
MaterialPropertyLib::VariableArray variables;
@@ -26,10 +25,23 @@ void PorosityModel::eval(SpaceTimeData const& x_t,
porosity_data.phi =
mpl_porosity.template value<double>(variables, x_t.x, x_t.t, x_t.dt);
+}
+
+void PorosityModel::dEval(SpaceTimeData const& x_t, MediaData const& media_data,
+ PorosityData const& porosity_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ PorosityDerivativeData& porosity_d_data) const
+{
+ MaterialPropertyLib::VariableArray variables;
+
+ auto const& mpl_porosity =
+ media_data.medium[MaterialPropertyLib::PropertyType::porosity];
porosity_d_data.dphi_dT = mpl_porosity.template dValue<double>(
variables, MaterialPropertyLib::Variable::temperature, x_t.x, x_t.t,
x_t.dt);
+
+ porosity_d_data.dphi_L_dp_cap = dS_L_dp_cap() * porosity_data.phi;
}
template <int DisplacementDim>
@@ -39,8 +51,7 @@ void PorosityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>::eval(
BiotData const& biot,
StrainData<DisplacementDim> const& strain_data,
SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
- PorosityData& porosity_data,
- PorosityDerivativeData& porosity_d_data) const
+ PorosityData& porosity_data) const
{
MaterialPropertyLib::VariableArray variables;
@@ -60,6 +71,31 @@ void PorosityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>::eval(
(1. + s_therm_exp_data.thermal_volume_strain - biot() * div_u);
porosity_data.phi = 1. - phi_S;
+}
+
+template <int DisplacementDim>
+void PorosityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>::dEval(
+ SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ PorosityData const& porosity_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ BiotData const& biot,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ StrainData<DisplacementDim> const& strain_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);
auto const dphi_0_dT = mpl_porosity.template dValue<double>(
variables, MaterialPropertyLib::Variable::temperature, x_t.x, x_t.t,
@@ -69,6 +105,8 @@ void PorosityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>::eval(
dphi_0_dT *
(1. + s_therm_exp_data.thermal_volume_strain - biot() * div_u) -
(1. - phi_0) * s_therm_exp_data.beta_T_SR;
+
+ porosity_d_data.dphi_L_dp_cap = dS_L_dp_cap() * porosity_data.phi;
}
template struct PorosityModelNonConstantSolidPhaseVolumeFraction<2>;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h
index f7e0d5d8f09a428005779f300cda2596ac9bd9f1..49884ed45eddd84dfe40308f5353bcda18538421 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h
@@ -13,6 +13,7 @@
#include "Biot.h"
#include "ProcessLib/ConstitutiveRelations/StrainData.h"
#include "ProcessLib/Reflection/ReflectionData.h"
+#include "Saturation.h"
#include "SolidThermalExpansion.h"
namespace ProcessLib::TH2M
@@ -22,6 +23,10 @@ namespace ConstitutiveRelations
struct PorosityDerivativeData
{
double dphi_dT = nan;
+ // dphi_G_dp_GR = -ds_L_dp_GR * phi = 0;
+ // dphi_L_dp_GR = ds_L_dp_GR * phi = 0;
+ double dphi_L_dp_cap = nan;
+ // dphi_G_dp_cap = -dphi_L_dp_cap
};
struct PorosityData
@@ -41,8 +46,13 @@ struct PorosityModel
{
void eval(SpaceTimeData const& x_t,
MediaData const& media_data,
- PorosityData& porosity_data,
- PorosityDerivativeData& porosity_d_data) const;
+ PorosityData& porosity_data) const;
+
+ void dEval(SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ PorosityData const& porosity_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ PorosityDerivativeData& porosity_d_data) const;
};
template <int DisplacementDim>
@@ -54,7 +64,16 @@ struct PorosityModelNonConstantSolidPhaseVolumeFraction
BiotData const& biot,
StrainData<DisplacementDim> const& strain_data,
SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
- PorosityData& porosity_data,
+ PorosityData& porosity_data) const;
+
+ void dEval(
+ SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ PorosityData const& porosity_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ BiotData const& biot,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ StrainData<DisplacementDim> const& strain_data,
PorosityDerivativeData& porosity_d_data) const;
};
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Saturation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Saturation.cpp
index 5f4693174a6e9b7f4cd0c88f2a5c9b3efbfd67ce..8ca73adaa1b7820caccfd10db4d4fb64f0fee18e 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Saturation.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Saturation.cpp
@@ -16,8 +16,7 @@ namespace ConstitutiveRelations
void SaturationModel::eval(SpaceTimeData const& x_t,
MediaData const& media_data,
CapillaryPressureData const& p_cap,
- SaturationData& S_L_data,
- SaturationDataDeriv& dS_L_data) const
+ SaturationData& S_L_data) const
{
namespace MPL = MaterialPropertyLib;
MPL::VariableArray variables;
@@ -27,6 +26,18 @@ void SaturationModel::eval(SpaceTimeData const& x_t,
S_L_data.S_L = medium.property(MPL::PropertyType::saturation)
.template value<double>(variables, x_t.x, x_t.t, x_t.dt);
+}
+
+void SaturationModel::dEval(SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ CapillaryPressureData const& p_cap,
+ SaturationDataDeriv& dS_L_data) const
+{
+ namespace MPL = MaterialPropertyLib;
+ MPL::VariableArray variables;
+ variables.capillary_pressure = p_cap();
+
+ auto const& medium = media_data.medium;
dS_L_data() = medium.property(MPL::PropertyType::saturation)
.template dValue<double>(
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Saturation.h b/ProcessLib/TH2M/ConstitutiveRelations/Saturation.h
index f8ab3071673e39bce24b1a67305bdc3551f8db17..66dbbea5dcc486d838923da8cc2130a3c8a13773 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Saturation.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Saturation.h
@@ -34,8 +34,12 @@ using SaturationDataDeriv =
struct SaturationModel
{
void eval(SpaceTimeData const& x_t, MediaData const& media_data,
- CapillaryPressureData const& p_cap, SaturationData& S_L_data,
- SaturationDataDeriv& dS_L_data) const;
+ CapillaryPressureData const& p_cap,
+ SaturationData& S_L_data) const;
+
+ void dEval(SpaceTimeData const& x_t, MediaData const& media_data,
+ CapillaryPressureData const& p_cap,
+ SaturationDataDeriv& dS_L_data) const;
};
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.cpp
index 9978c1516908253b798ee355920cc70fb4c5809b..5a764bc2015c29aef9c397c166ae62042a9eae38 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.cpp
@@ -14,11 +14,25 @@ namespace ProcessLib::TH2M
namespace ConstitutiveRelations
{
-void SolidDensityModel::eval(
+void SolidDensityModel::eval(SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ TemperatureData const& T_data,
+ SolidDensityData& solid_density_data) const
+{
+ MaterialPropertyLib::VariableArray variables;
+ variables.temperature = T_data.T;
+
+ auto const& mpl_solid_density =
+ media_data.solid[MaterialPropertyLib::PropertyType::density];
+
+ solid_density_data.rho_SR = mpl_solid_density.template value<double>(
+ variables, x_t.x, x_t.t, x_t.dt);
+}
+
+void SolidDensityModel::dEval(
SpaceTimeData const& x_t,
MediaData const& media_data,
TemperatureData const& T_data,
- SolidDensityData& solid_density_data,
SolidDensityDerivativeData& solid_density_d_data) const
{
MaterialPropertyLib::VariableArray variables;
@@ -27,9 +41,6 @@ void SolidDensityModel::eval(
auto const& mpl_solid_density =
media_data.solid[MaterialPropertyLib::PropertyType::density];
- solid_density_data.rho_SR = mpl_solid_density.template value<double>(
- variables, x_t.x, x_t.t, x_t.dt);
-
solid_density_d_data.drho_SR_dT = mpl_solid_density.template dValue<double>(
variables, MaterialPropertyLib::Variable::temperature, x_t.x, x_t.t,
x_t.dt);
@@ -43,8 +54,7 @@ void SolidDensityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>::
BiotData const& biot,
StrainData<DisplacementDim> const& strain_data,
SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
- SolidDensityData& solid_density_data,
- SolidDensityDerivativeData& solid_density_d_data) const
+ SolidDensityData& solid_density_data) const
{
MaterialPropertyLib::VariableArray variables;
variables.temperature = T_data.T;
@@ -63,6 +73,31 @@ void SolidDensityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>::
solid_density_data.rho_SR =
rho_ref_SR *
(1. - s_therm_exp_data.thermal_volume_strain + (biot() - 1.) * div_u);
+}
+
+template <int DisplacementDim>
+void SolidDensityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>::
+ dEval(SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ TemperatureData const& T_data,
+ BiotData const& biot,
+ StrainData<DisplacementDim> const& strain_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ SolidDensityDerivativeData& solid_density_d_data) const
+{
+ MaterialPropertyLib::VariableArray variables;
+ variables.temperature = T_data.T;
+
+ 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);
+
+ auto const& mpl_solid_density =
+ media_data.solid[MaterialPropertyLib::PropertyType::density];
+
+ auto const rho_ref_SR = mpl_solid_density.template value<double>(
+ variables, x_t.x, x_t.t, x_t.dt);
solid_density_d_data.drho_SR_dT =
mpl_solid_density.template dValue<double>(
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h
index ad14ab20ff85c3de00210a6740324b59e3866e4f..d9f2d4890c4766c1dc1390db965a4b159f2d71d9 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h
@@ -43,8 +43,12 @@ struct SolidDensityModel
void eval(SpaceTimeData const& x_t,
MediaData const& media_data,
TemperatureData const& T_data,
- SolidDensityData& solid_density_data,
- SolidDensityDerivativeData& solid_density_d_data) const;
+ SolidDensityData& solid_density_data) const;
+
+ void dEval(SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ TemperatureData const& T_data,
+ SolidDensityDerivativeData& solid_density_d_data) const;
};
template <int DisplacementDim>
@@ -57,7 +61,15 @@ struct SolidDensityModelNonConstantSolidPhaseVolumeFraction
BiotData const& biot,
StrainData<DisplacementDim> const& strain_data,
SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
- SolidDensityData& solid_density_data,
+ SolidDensityData& solid_density_data) const;
+
+ void dEval(
+ SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ TemperatureData const& T_data,
+ BiotData const& biot,
+ StrainData<DisplacementDim> const& strain_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
SolidDensityDerivativeData& solid_density_d_data) const;
};
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/TEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/TEquation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0dd5b7b814ba1b5b05c67214560866c3c35904bc
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/TEquation.cpp
@@ -0,0 +1,147 @@
+/**
+ * \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 "TEquation.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+void FT1Model::eval(
+ double const dt,
+ InternalEnergyData const& internal_energy_data,
+ PrevState<InternalEnergyData> const& internal_energy_data_prev,
+ FT1Data& fT_1) const
+{
+ if (dt == 0)
+ {
+ fT_1.m = 0;
+ return;
+ }
+
+ auto const rho_u_eff_dot =
+ (internal_energy_data() - **internal_energy_data_prev) / dt;
+ fT_1.m = rho_u_eff_dot;
+}
+
+void FT1Model::dEval(double const dt,
+ EffectiveVolumetricInternalEnergyDerivatives const&
+ effective_volumetric_internal_energy_d_data,
+ FT1DerivativeData& dfT_1) const
+{
+ if (dt == 0)
+ {
+ dfT_1.dp_GR = 0;
+ dfT_1.dp_cap = 0;
+ dfT_1.dT = 0;
+ return;
+ }
+
+ dfT_1.dp_GR =
+ effective_volumetric_internal_energy_d_data.drho_u_eff_dp_GR / dt;
+
+ dfT_1.dp_cap =
+ effective_volumetric_internal_energy_d_data.drho_u_eff_dp_cap / dt;
+
+ dfT_1.dT = effective_volumetric_internal_energy_d_data.drho_u_eff_dT / dt;
+}
+
+template <int DisplacementDim>
+void FT2Model<DisplacementDim>::eval(
+ DarcyVelocityData<DisplacementDim> const& darcy_velocity_data,
+ FluidDensityData const& fluid_density_data,
+ FluidEnthalpyData const& fluid_enthalpy_data,
+ FT2Data<DisplacementDim>& fT_2) const
+{
+ fT_2.A.noalias() = fluid_density_data.rho_GR * fluid_enthalpy_data.h_G *
+ darcy_velocity_data.w_GS +
+ fluid_density_data.rho_LR * fluid_enthalpy_data.h_L *
+ darcy_velocity_data.w_LS;
+}
+
+template <int DisplacementDim>
+void FT2Model<DisplacementDim>::dEval(
+ DarcyVelocityData<DisplacementDim> const& darcy_velocity_data,
+ FluidDensityData const& fluid_density_data,
+ FluidEnthalpyData const& fluid_enthalpy_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ PhaseTransitionData const& phase_transition_data,
+ SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+ ViscosityData const& viscosity_data,
+ FT2DerivativeData<DisplacementDim>& dfT_2) const
+{
+ auto const k_over_mu_G =
+ permeability_data.Ki * permeability_data.k_rel_G / viscosity_data.mu_GR;
+ auto const k_over_mu_L =
+ permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR;
+
+ dfT_2.dp_GR_Npart = phase_transition_data.drho_GR_dp_GR *
+ fluid_enthalpy_data.h_G * darcy_velocity_data.w_GS +
+ fluid_density_data.rho_GR * fluid_enthalpy_data.h_G *
+ k_over_mu_G * phase_transition_data.drho_GR_dp_GR *
+ specific_body_force();
+ dfT_2.dp_GR_gradNpart =
+ fluid_density_data.rho_GR * fluid_enthalpy_data.h_G * k_over_mu_G -
+ fluid_density_data.rho_LR * fluid_enthalpy_data.h_L * k_over_mu_L;
+
+ // From p_LR = p_GR - p_cap it follows for
+ // drho_LR/dp_GR = drho_LR/dp_LR * dp_LR/dp_GR
+ // = drho_LR/dp_LR * (dp_GR/dp_GR - dp_cap/dp_GR)
+ // = drho_LR/dp_LR * (1 - 0)
+ double const drho_LR_dp_cap = -phase_transition_data.drho_LR_dp_LR;
+
+ dfT_2.dp_cap_Npart =
+ -drho_LR_dp_cap * fluid_enthalpy_data.h_L * darcy_velocity_data.w_LS -
+ fluid_density_data.rho_LR * fluid_enthalpy_data.h_L * k_over_mu_L *
+ drho_LR_dp_cap * specific_body_force();
+ dfT_2.dp_cap_gradNpart =
+ fluid_density_data.rho_LR * fluid_enthalpy_data.h_L * k_over_mu_L;
+
+ dfT_2.dT = phase_transition_data.drho_GR_dT * fluid_enthalpy_data.h_G *
+ darcy_velocity_data.w_GS +
+ fluid_density_data.rho_GR * phase_transition_data.dh_G_dT *
+ darcy_velocity_data.w_GS +
+ phase_transition_data.drho_LR_dT * fluid_enthalpy_data.h_L *
+ darcy_velocity_data.w_LS +
+ fluid_density_data.rho_LR * phase_transition_data.dh_L_dT *
+ darcy_velocity_data.w_LS;
+ // TODO (naumov) + k_over_mu_G * drho_GR_dT * b + k_over_mu_L *
+ // drho_LR_dT * b
+}
+
+template struct FT2Model<2>;
+template struct FT2Model<3>;
+
+template <int DisplacementDim>
+void FT3Model<DisplacementDim>::eval(
+ ConstituentDensityData const& constituent_density_data,
+ DarcyVelocityData<DisplacementDim> const& darcy_velocity_data,
+ DiffusionVelocityData<DisplacementDim> const& diffusion_velocity_data,
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+ FT3Data<DisplacementDim>& fT_3) const
+{
+ fT_3.N =
+ (fluid_density_data.rho_GR * darcy_velocity_data.w_GS.transpose() +
+ fluid_density_data.rho_LR * darcy_velocity_data.w_LS.transpose()) *
+ specific_body_force();
+
+ fT_3.gradN.noalias() =
+ constituent_density_data.rho_C_GR * phase_transition_data.hCG *
+ diffusion_velocity_data.d_CG +
+ constituent_density_data.rho_W_GR * phase_transition_data.hWG *
+ diffusion_velocity_data.d_WG;
+}
+
+template struct FT3Model<2>;
+template struct FT3Model<3>;
+
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/TEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/TEquation.h
new file mode 100644
index 0000000000000000000000000000000000000000..94d80aaf05f43f1ae70209aaf1862be10edd4d21
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/TEquation.h
@@ -0,0 +1,112 @@
+/**
+ * \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"
+#include "ConstitutiveDensity.h"
+#include "DarcyVelocity.h"
+#include "DiffusionVelocity.h"
+#include "Enthalpy.h"
+#include "FluidDensity.h"
+#include "InternalEnergy.h"
+#include "PermeabilityData.h"
+#include "Viscosity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+struct FT1Data
+{
+ double m = nan;
+};
+
+struct FT1DerivativeData
+{
+ double dp_GR = nan;
+ double dp_cap = nan;
+ double dT = nan;
+};
+
+struct FT1Model
+{
+ void eval(double const dt,
+ InternalEnergyData const& internal_energy_data,
+ PrevState<InternalEnergyData> const& internal_energy_data_prev,
+ FT1Data& fT_1) const;
+
+ void dEval(double const dt,
+ EffectiveVolumetricInternalEnergyDerivatives const&
+ effective_volumetric_internal_energy_d_data,
+ FT1DerivativeData& dfT_1) const;
+};
+
+template <int DisplacementDim>
+struct FT2Data
+{
+ GlobalDimVector<DisplacementDim> A;
+};
+
+template <int DisplacementDim>
+struct FT2DerivativeData
+{
+ GlobalDimVector<DisplacementDim> dp_GR_Npart;
+ GlobalDimMatrix<DisplacementDim> dp_GR_gradNpart;
+ GlobalDimVector<DisplacementDim> dp_cap_Npart;
+ GlobalDimMatrix<DisplacementDim> dp_cap_gradNpart;
+ GlobalDimVector<DisplacementDim> dT;
+};
+
+template <int DisplacementDim>
+struct FT2Model
+{
+ void eval(DarcyVelocityData<DisplacementDim> const& darcy_velocity_data,
+ FluidDensityData const& fluid_density_data,
+ FluidEnthalpyData const& fluid_enthalpy_data,
+ FT2Data<DisplacementDim>& fT_2) const;
+
+ void dEval(
+ DarcyVelocityData<DisplacementDim> const& darcy_velocity_data,
+ FluidDensityData const& fluid_density_data,
+ FluidEnthalpyData const& fluid_enthalpy_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ PhaseTransitionData const& phase_transition_data,
+ SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+ ViscosityData const& viscosity_data,
+ FT2DerivativeData<DisplacementDim>& dfT_2) const;
+};
+
+extern template struct FT2Model<2>;
+extern template struct FT2Model<3>;
+
+template <int DisplacementDim>
+struct FT3Data
+{
+ double N = nan;
+ GlobalDimVector<DisplacementDim> gradN;
+};
+
+template <int DisplacementDim>
+struct FT3Model
+{
+ void eval(
+ ConstituentDensityData const& constituent_density_data,
+ DarcyVelocityData<DisplacementDim> const& darcy_velocity_data,
+ DiffusionVelocityData<DisplacementDim> const& diffusion_velocity_data,
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+ FT3Data<DisplacementDim>& fT_3) const;
+};
+
+extern template struct FT3Model<2>;
+extern template struct FT3Model<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp
index abd9e0cc9f3f4d2bdf7cd19f3627209e21b28edd..358ee797cc7005035ee0fdeeb6a1df98cca2f271 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp
@@ -19,8 +19,7 @@ template <int DisplacementDim>
void ThermalConductivityModel<DisplacementDim>::eval(
SpaceTimeData const& x_t, MediaData const& media_data,
TemperatureData const& T_data, PorosityData const& porosity_data,
- PorosityDerivativeData const& porosity_d_data,
- SaturationData const& S_L_data, SaturationDataDeriv const& dS_L_dp_cap,
+ SaturationData const& S_L_data,
ThermalConductivityData<DisplacementDim>& thermal_conductivity_data) const
{
namespace MPL = MaterialPropertyLib;
@@ -34,6 +33,22 @@ void ThermalConductivityModel<DisplacementDim>::eval(
thermal_conductivity_data.lambda = MPL::formEigenTensor<DisplacementDim>(
mpl_thermal_conductivity.value(variables, x_t.x, x_t.t, x_t.dt));
+}
+
+template <int DisplacementDim>
+void ThermalConductivityModel<DisplacementDim>::dEval(
+ SpaceTimeData const& x_t, MediaData const& media_data,
+ TemperatureData const& T_data, PorosityData const& porosity_data,
+ PorosityDerivativeData const& porosity_d_data,
+ SaturationData const& S_L_data,
+ ThermalConductivityDerivativeData<DisplacementDim>&
+ thermal_conductivity_d_data) const
+{
+ namespace MPL = MaterialPropertyLib;
+ MPL::VariableArray variables;
+ variables.temperature = T_data.T;
+ variables.porosity = porosity_data.phi;
+ variables.liquid_saturation = S_L_data.S_L;
// Derivatives computed here and not in the MPL property because various
// derivatives are not available in the VariableArray.
@@ -83,13 +98,9 @@ void ThermalConductivityModel<DisplacementDim>::eval(
x_t.t, x_t.dt))
: MPL::formEigenTensor<DisplacementDim>(0.);
- // dphi_G_dp_GR = -ds_L_dp_GR * phi = 0;
- double const dphi_G_dp_cap = -dS_L_dp_cap() * porosity_data.phi;
- // dphi_L_dp_GR = ds_L_dp_GR * phi = 0;
- double const dphi_L_dp_cap = dS_L_dp_cap() * porosity_data.phi;
-
- thermal_conductivity_data.dlambda_dp_cap =
- dphi_G_dp_cap * lambdaGR + dphi_L_dp_cap * lambdaLR;
+ thermal_conductivity_d_data.dlambda_dp_cap =
+ -porosity_d_data.dphi_L_dp_cap * lambdaGR +
+ porosity_d_data.dphi_L_dp_cap * lambdaLR;
double const phi_L = S_L_data.S_L * porosity_data.phi;
double const phi_G = (1. - S_L_data.S_L) * porosity_data.phi;
@@ -99,13 +110,14 @@ void ThermalConductivityModel<DisplacementDim>::eval(
// dphi_G_dT = -dS_L/dT * phi + (1 - S_L) * dphi_dT = (1 - S_L) * dphi_dT
// dphi_L_dT = dS_L/dT * phi + S_L * dphi_dT S_L * dphi_dT
// dphi_S_dT = -dphi_dT -dphi_dT
- thermal_conductivity_data.dlambda_dT =
+ thermal_conductivity_d_data.dlambda_dT =
(1 - S_L_data.S_L) * porosity_d_data.dphi_dT * lambdaGR +
phi_G * dlambda_GR_dT +
S_L_data.S_L * porosity_d_data.dphi_dT * lambdaLR +
+phi_L * dlambda_LR_dT - porosity_d_data.dphi_dT * lambdaSR +
phi_S * dlambda_SR_dT;
}
+
template struct ThermalConductivityModel<2>;
template struct ThermalConductivityModel<3>;
} // namespace ConstitutiveRelations
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h
index fbf34021e55633186144f6fcc267acf6e5d1f45c..35720adac3fee11f0e318e8111151868261f8cc8 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h
@@ -22,6 +22,11 @@ template <int DisplacementDim>
struct ThermalConductivityData
{
GlobalDimMatrix<DisplacementDim> lambda;
+};
+
+template <int DisplacementDim>
+struct ThermalConductivityDerivativeData
+{
// Currently unused, but there is a comment in TH2MFEM-impl.h referring to
// this matrix
// GlobalDimMatrix<DisplacementDim> dlambda_dp_GR;
@@ -34,11 +39,16 @@ struct ThermalConductivityModel
{
void eval(SpaceTimeData const& x_t, MediaData const& media_data,
TemperatureData const& T_data, PorosityData const& porosity_data,
- PorosityDerivativeData const& porosity_d_data,
SaturationData const& S_L_data,
- SaturationDataDeriv const& dS_L_dp_cap,
ThermalConductivityData<DisplacementDim>&
thermal_conductivity_data) const;
+
+ void dEval(SpaceTimeData const& x_t, MediaData const& media_data,
+ TemperatureData const& T_data, PorosityData const& porosity_data,
+ PorosityDerivativeData const& porosity_d_data,
+ SaturationData const& S_L_data,
+ ThermalConductivityDerivativeData<DisplacementDim>&
+ thermal_conductivity_d_data) const;
};
extern template struct ThermalConductivityModel<2>;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dcea844e0a42ee7603f51c5099bcde21c6170cab
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp
@@ -0,0 +1,47 @@
+/**
+ * \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 "UEquation.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void FU1KUTModel<DisplacementDim>::dEval(
+ SolidMechanicsDataStateless<DisplacementDim> const& s_mech_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ FU1KUTDerivativeData<DisplacementDim>& dfu_1_KuT) const
+{
+ dfu_1_KuT.dT = s_mech_data.stiffness_tensor *
+ s_therm_exp_data.solid_linear_thermal_expansivity_vector;
+}
+
+template struct FU1KUTModel<2>;
+template struct FU1KUTModel<3>;
+
+void FU2KUpCModel::eval(BiotData const& biot_data,
+ BishopsData const& chi_S_L,
+ FU2KUpCData& fu_2_KupC) const
+{
+ fu_2_KupC.m = biot_data() * chi_S_L.chi_S_L;
+}
+
+void FU2KUpCModel::dEval(BiotData const& biot_data,
+ BishopsData const& chi_S_L,
+ CapillaryPressureData const& p_cap,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ FU2KUpCDerivativeData& dfu_2_KupC) const
+{
+ dfu_2_KupC.dp_cap =
+ biot_data() * chi_S_L.dchi_dS_L * dS_L_dp_cap() * p_cap();
+}
+
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h
new file mode 100644
index 0000000000000000000000000000000000000000..70388bf50d2228018e76958c5aa2898b45dbdaa2
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h
@@ -0,0 +1,64 @@
+/**
+ * \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"
+#include "Biot.h"
+#include "Bishops.h"
+#include "Saturation.h"
+#include "SolidMechanics.h"
+#include "SolidThermalExpansion.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+struct FU1KUTDerivativeData
+{
+ KelvinVector<DisplacementDim> dT;
+};
+
+template <int DisplacementDim>
+struct FU1KUTModel
+{
+ void dEval(
+ SolidMechanicsDataStateless<DisplacementDim> const& s_mech_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ FU1KUTDerivativeData<DisplacementDim>& dfu_1_KuT) const;
+};
+
+extern template struct FU1KUTModel<2>;
+extern template struct FU1KUTModel<3>;
+
+struct FU2KUpCData
+{
+ double m = nan;
+};
+
+struct FU2KUpCDerivativeData
+{
+ double dp_cap = nan;
+};
+
+struct FU2KUpCModel
+{
+ void eval(BiotData const& biot_data,
+ BishopsData const& chi_S_L,
+ FU2KUpCData& fu_2_KupC) const;
+
+ void dEval(BiotData const& biot_data,
+ BishopsData const& chi_S_L,
+ CapillaryPressureData const& p_cap,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ FU2KUpCDerivativeData& dfu_2_KupC) const;
+};
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..370e22cd5d0430c60a550a03b1de1ea80d1d40a6
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
@@ -0,0 +1,449 @@
+/**
+ * \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 "WEquation.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void FW1Model<DisplacementDim>::eval(
+ AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ FW1Data<DisplacementDim>& fW_1) const
+{
+ fW_1.A = advection_data.advection_W_G * fluid_density_data.rho_GR +
+ advection_data.advection_W_L * fluid_density_data.rho_LR;
+}
+
+template struct FW1Model<2>;
+template struct FW1Model<3>;
+
+void FW2Model::eval(BiotData const biot_data,
+ CapillaryPressureData const pCap,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const beta_p_SR,
+ FW2Data& fW_2) const
+{
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ double const rho_W_FR =
+ S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR();
+
+ fW_2.a =
+ porosity_data.phi * (rho_W_LR() - constituent_density_data.rho_W_GR) -
+ rho_W_FR * pCap() * (biot_data() - porosity_data.phi) * beta_p_SR();
+}
+
+void FW2Model::dEval(BiotData const& biot_data,
+ CapillaryPressureData const pCap,
+ ConstituentDensityData const& constituent_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PorosityDerivativeData const& porosity_d_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ SolidCompressibilityData const& beta_p_SR,
+ FW2DerivativeData& dfW_2) const
+{
+ double const S_L = S_L_data.S_L;
+ double const S_G = 1. - S_L;
+
+ double const drho_C_FR_dp_GR =
+ /*(dS_G_dp_GR = 0) * constituent_density_data.rho_C_GR +*/
+ S_G * phase_transition_data.drho_C_GR_dp_GR +
+ /*(dS_L_dp_GR = 0) * constituent_density_data.rho_C_LR +*/
+ S_L * phase_transition_data.drho_C_LR_dp_GR;
+
+ dfW_2.dp_GR = -porosity_data.phi * phase_transition_data.drho_C_GR_dp_GR -
+ drho_C_FR_dp_GR * pCap() * (biot_data() - porosity_data.phi) *
+ beta_p_SR();
+
+ double const dfW_2a_dp_cap =
+ porosity_data.phi * (-phase_transition_data.drho_W_LR_dp_LR -
+ phase_transition_data.drho_W_GR_dp_cap);
+ double const rho_W_FR =
+ S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR();
+
+ double const drho_W_FR_dp_cap =
+ -dS_L_dp_cap() * constituent_density_data.rho_W_GR +
+ S_G * phase_transition_data.drho_W_GR_dp_cap +
+ dS_L_dp_cap() * rho_W_LR() -
+ S_L * phase_transition_data.drho_W_LR_dp_LR;
+
+ double const dfW_2b_dp_cap =
+ drho_W_FR_dp_cap * pCap() * (biot_data() - porosity_data.phi) *
+ beta_p_SR() +
+ rho_W_FR * (biot_data() - porosity_data.phi) * beta_p_SR();
+
+ dfW_2.dp_cap = dfW_2a_dp_cap - dfW_2b_dp_cap;
+
+ double const drho_W_FR_dT = S_G * phase_transition_data.drho_W_GR_dT +
+ S_L * phase_transition_data.drho_W_LR_dT;
+
+ double const dfW_2a_dT =
+ porosity_d_data.dphi_dT *
+ (rho_W_LR() - constituent_density_data.rho_W_GR) +
+ porosity_data.phi * (phase_transition_data.drho_W_LR_dT -
+ phase_transition_data.drho_W_GR_dT);
+ double const dfW_2b_dT =
+ drho_W_FR_dT * pCap() * (biot_data() - porosity_data.phi) *
+ beta_p_SR() -
+ rho_W_FR * pCap() * porosity_d_data.dphi_dT * beta_p_SR();
+
+ dfW_2.dT = dfW_2a_dT - dfW_2b_dT;
+}
+
+void FW3aModel::eval(
+ double const dt,
+ ConstituentDensityData const& constituent_density_data,
+ PrevState<ConstituentDensityData> const& constituent_density_data_prev,
+ PrevState<PureLiquidDensityData> const& rho_W_LR_prev,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ FW3aData& fW_3a) const
+{
+ if (dt == 0.)
+ {
+ fW_3a.a = 0;
+ return;
+ }
+
+ double const rho_W_GR_dot = (constituent_density_data.rho_W_GR -
+ constituent_density_data_prev->rho_W_GR) /
+ dt;
+ double const rho_W_LR_dot = (rho_W_LR() - **rho_W_LR_prev) / dt;
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ fW_3a.a = S_G * rho_W_GR_dot + S_L * rho_W_LR_dot;
+}
+
+void FW3aModel::dEval(
+ double const dt,
+ ConstituentDensityData const& constituent_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PrevState<ConstituentDensityData> const& constituent_density_data_prev,
+ PrevState<PureLiquidDensityData> const& rho_W_LR_prev,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ FW3aDerivativeData& dfW_3a) const
+{
+ if (dt == 0.)
+ {
+ dfW_3a.dp_GR = 0.;
+ dfW_3a.dp_cap = 0.;
+ dfW_3a.dT = 0.;
+ return;
+ }
+
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+
+ double const rho_W_GR_dot = (constituent_density_data.rho_W_GR -
+ constituent_density_data_prev->rho_W_GR) /
+ dt;
+ double const rho_W_LR_dot = (rho_W_LR() - **rho_W_LR_prev) / dt;
+
+ dfW_3a.dp_GR = /*(ds_G_dp_GR = 0) * rho_W_GR_dot +*/
+ S_G * phase_transition_data.drho_W_GR_dp_GR / dt +
+ /*(ds_L_dp_GR = 0) * rho_W_LR_dot +*/
+ S_L * phase_transition_data.drho_W_LR_dp_GR / dt;
+
+ dfW_3a.dp_cap = -dS_L_dp_cap() * rho_W_GR_dot +
+ S_G * phase_transition_data.drho_W_GR_dp_cap / dt +
+ dS_L_dp_cap() * rho_W_LR_dot -
+ S_L * phase_transition_data.drho_W_LR_dp_LR / dt;
+ dfW_3a.dT = S_G * phase_transition_data.drho_W_GR_dT / dt +
+ S_L * phase_transition_data.drho_W_LR_dT / dt;
+}
+
+template <int DisplacementDim>
+void FW4LWpGModel<DisplacementDim>::eval(
+ AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FW4LWpGData<DisplacementDim>& fW_4_LWpG) const
+{
+ GlobalDimMatrix<DisplacementDim> const advection_W =
+ advection_data.advection_W_G + advection_data.advection_W_L;
+
+ double const sD_G = phase_transition_data.diffusion_coefficient_vapour;
+ double const sD_L = phase_transition_data.diffusion_coefficient_solute;
+
+ double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi;
+ double const phi_L = S_L_data.S_L * porosity_data.phi;
+
+ double const diffusion_WGpGR = phi_G * fluid_density_data.rho_GR * sD_G *
+ phase_transition_data.dxmWG_dpGR;
+ double const diffusion_WLpGR = phi_L * fluid_density_data.rho_LR * sD_L *
+ phase_transition_data.dxmWL_dpGR;
+ double const diffusion_W_pGR = diffusion_WGpGR + diffusion_WLpGR;
+
+ auto const I =
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
+ fW_4_LWpG.L.noalias() = diffusion_W_pGR * I + advection_W;
+}
+
+template <int DisplacementDim>
+void FW4LWpGModel<DisplacementDim>::dEval(
+ ConstituentDensityData const& constituent_density_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ PhaseTransitionData const& phase_transition_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ ViscosityData const& viscosity_data,
+ FW4LWpGDerivativeData<DisplacementDim>& dfW_4_LWpG) const
+{
+ ////// Diffusion Part /////
+ // TODO (naumov) d(diffusion_W_p)/dX for dxmW*/d* != 0
+
+ auto const k_over_mu_G =
+ permeability_data.Ki * permeability_data.k_rel_G / viscosity_data.mu_GR;
+ auto const k_over_mu_L =
+ permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR;
+
+ // dk_over_mu_G_dp_GR = ip_out.permeability_data.Ki *
+ // ip_out.permeability_data.dk_rel_G_dS_L *
+ // (ds_L_dp_GR = 0) /
+ // ip_cv.viscosity_data.mu_GR = 0;
+ // dk_over_mu_L_dp_GR = ip_out.permeability_data.Ki *
+ // ip_out.permeability_data.dk_rel_L_dS_L *
+ // (ds_L_dp_GR = 0) /
+ // ip_cv.viscosity_data.mu_LR = 0;
+ auto const dk_over_mu_G_dp_cap = permeability_data.Ki *
+ permeability_data.dk_rel_G_dS_L *
+ dS_L_dp_cap() / viscosity_data.mu_GR;
+
+ auto const dk_over_mu_L_dp_cap = permeability_data.Ki *
+ permeability_data.dk_rel_L_dS_L *
+ dS_L_dp_cap() / viscosity_data.mu_LR;
+
+ dfW_4_LWpG.dp_GR = phase_transition_data.drho_W_GR_dp_GR * k_over_mu_G
+ // + rhoWGR * (dk_over_mu_G_dp_GR = 0)
+ + phase_transition_data.drho_W_LR_dp_GR * k_over_mu_L
+ // + rhoWLR * (dk_over_mu_L_dp_GR = 0)
+ ;
+
+ dfW_4_LWpG.dp_cap =
+ phase_transition_data.drho_W_GR_dp_cap * k_over_mu_G +
+ constituent_density_data.rho_W_GR * dk_over_mu_G_dp_cap +
+ -phase_transition_data.drho_W_LR_dp_LR * k_over_mu_L +
+ rho_W_LR() * dk_over_mu_L_dp_cap;
+
+ dfW_4_LWpG.dT = phase_transition_data.drho_W_GR_dT * k_over_mu_G
+ //+ rhoWGR * (dk_over_mu_G_dT != 0 TODO for mu_G(T))
+ + phase_transition_data.drho_W_LR_dT * k_over_mu_L
+ //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_G(T))
+ ;
+}
+
+template struct FW4LWpGModel<2>;
+template struct FW4LWpGModel<3>;
+
+template <int DisplacementDim>
+void FW4LWpCModel<DisplacementDim>::eval(
+ AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FW4LWpCData<DisplacementDim>& fW_4_LWpC) const
+{
+ double const sD_G = phase_transition_data.diffusion_coefficient_vapour;
+ double const sD_L = phase_transition_data.diffusion_coefficient_solute;
+
+ double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi;
+ double const phi_L = S_L_data.S_L * porosity_data.phi;
+
+ double const diffusion_WGpCap = phi_G * fluid_density_data.rho_GR * sD_G *
+ phase_transition_data.dxmWG_dpCap;
+ double const diffusion_WLpCap = phi_L * fluid_density_data.rho_LR * sD_L *
+ phase_transition_data.dxmWL_dpCap;
+
+ double const diffusion_W_pCap = diffusion_WGpCap + diffusion_WLpCap;
+
+ auto const I =
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
+
+ fW_4_LWpC.L.noalias() = diffusion_W_pCap * I - advection_data.advection_W_L;
+}
+
+template <int DisplacementDim>
+void FW4LWpCModel<DisplacementDim>::dEval(
+ AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ ViscosityData const& viscosity_data,
+ FW4LWpCDerivativeData<DisplacementDim>& dfW_4_LWpC) const
+{
+ ////// Diffusion Part /////
+ // TODO (naumov) d(diffusion_W_pCap)/dX for dxmW*/d* != 0
+
+ ////// Advection part /////
+ GlobalDimMatrix<DisplacementDim> const k_over_mu_L =
+ permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR;
+
+ dfW_4_LWpC.dp_GR = phase_transition_data.drho_W_LR_dp_GR * k_over_mu_L
+ //+ rhoWLR * (dk_over_mu_L_dp_GR = 0)
+ ;
+
+ double const sD_G = phase_transition_data.diffusion_coefficient_vapour;
+ double const sD_L = phase_transition_data.diffusion_coefficient_solute;
+
+ double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi;
+ double const phi_L = S_L_data.S_L * porosity_data.phi;
+
+ double const diffusion_WGpCap = phi_G * fluid_density_data.rho_GR * sD_G *
+ phase_transition_data.dxmWG_dpCap;
+ double const diffusion_WLpCap = phi_L * fluid_density_data.rho_LR * sD_L *
+ phase_transition_data.dxmWL_dpCap;
+
+ double const diffusion_W_pCap = diffusion_WGpCap + diffusion_WLpCap;
+
+ auto const I =
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
+
+ dfW_4_LWpC.dp_cap = diffusion_W_pCap * I - advection_data.advection_W_L;
+
+ auto const dk_over_mu_L_dp_cap = permeability_data.Ki *
+ permeability_data.dk_rel_L_dS_L *
+ dS_L_dp_cap() / viscosity_data.mu_LR;
+ dfW_4_LWpC.dp_cap = -phase_transition_data.drho_W_LR_dp_LR * k_over_mu_L +
+ rho_W_LR() * dk_over_mu_L_dp_cap;
+
+ dfW_4_LWpC.dT = phase_transition_data.drho_W_LR_dT * k_over_mu_L
+ //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T))
+ ;
+}
+
+template struct FW4LWpCModel<2>;
+template struct FW4LWpCModel<3>;
+
+template <int DisplacementDim>
+void FW4LWTModel<DisplacementDim>::eval(
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FW4LWTData<DisplacementDim>& fW_4_LWT) const
+{
+ double const sD_G = phase_transition_data.diffusion_coefficient_vapour;
+ double const sD_L = phase_transition_data.diffusion_coefficient_solute;
+
+ double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi;
+ double const phi_L = S_L_data.S_L * porosity_data.phi;
+
+ double const diffusion_W_G_T = phi_G * fluid_density_data.rho_GR * sD_G *
+ phase_transition_data.dxmWG_dT;
+ double const diffusion_W_L_T = phi_L * fluid_density_data.rho_LR * sD_L *
+ phase_transition_data.dxmWL_dT;
+
+ double const diffusion_W_T = diffusion_W_G_T + diffusion_W_L_T;
+
+ auto const I =
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
+
+ fW_4_LWT.L.noalias() = diffusion_W_T * I;
+}
+
+template struct FW4LWTModel<2>;
+template struct FW4LWTModel<3>;
+
+void FW4MWpGModel::eval(BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const& beta_p_SR,
+ FW4MWpGData& fW_4_MWpG) const
+{
+ double const S_L = S_L_data.S_L;
+ double const S_G = 1. - S_L;
+
+ double const rho_W_FR =
+ S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR();
+
+ fW_4_MWpG.m = rho_W_FR * (biot_data() - porosity_data.phi) * beta_p_SR();
+}
+
+void FW4MWpCModel::eval(BiotData const& biot_data,
+ CapillaryPressureData const pCap,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ PrevState<SaturationData> const& S_L_data_prev,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const& beta_p_SR,
+ FW4MWpCData& fW_4_MWpC) const
+{
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ double const rho_W_FR =
+ S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR();
+
+ fW_4_MWpC.m =
+ -rho_W_FR * (biot_data() - porosity_data.phi) * beta_p_SR() * S_L;
+
+ fW_4_MWpC.ml =
+ (porosity_data.phi * (rho_W_LR() - constituent_density_data.rho_W_GR) -
+ rho_W_FR * pCap() * (biot_data() - porosity_data.phi) * beta_p_SR()) *
+ (S_L - S_L_data_prev->S_L);
+}
+
+template <int DisplacementDim>
+void FW4MWTModel<DisplacementDim>::eval(
+ BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ FW4MWTData& fW_4_MWT) const
+{
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ double const rho_W_FR =
+ S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR();
+
+ fW_4_MWT.m = -rho_W_FR * (biot_data() - porosity_data.phi) *
+ s_therm_exp_data.beta_T_SR;
+}
+
+template struct FW4MWTModel<2>;
+template struct FW4MWTModel<3>;
+
+void FW4MWuModel::eval(BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ FW4MWuData& fW_4_MWu) const
+{
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ double const rho_W_FR =
+ S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR();
+
+ fW_4_MWu.m = rho_W_FR * biot_data();
+}
+
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
new file mode 100644
index 0000000000000000000000000000000000000000..09f1b09fc444375f154a8860589b1b37065f46f1
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
@@ -0,0 +1,277 @@
+/**
+ * \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 "Advection.h"
+#include "Base.h"
+#include "Biot.h"
+#include "ConstitutiveDensity.h"
+#include "FluidDensity.h"
+#include "PermeabilityData.h"
+#include "PhaseTransitionData.h"
+#include "Porosity.h"
+#include "Saturation.h"
+#include "SolidCompressibility.h"
+#include "Viscosity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+struct FW1Data
+{
+ GlobalDimMatrix<DisplacementDim> A;
+};
+
+template <int DisplacementDim>
+struct FW1Model
+{
+ void eval(AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ FW1Data<DisplacementDim>& fW_1) const;
+};
+
+extern template struct FW1Model<2>;
+extern template struct FW1Model<3>;
+
+struct FW2Data
+{
+ double a = nan;
+};
+
+struct FW2DerivativeData
+{
+ double dp_GR = nan;
+ double dp_cap = nan;
+ double dT = nan;
+};
+
+struct FW2Model
+{
+ void eval(BiotData const biot_data,
+ CapillaryPressureData const pCap,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const beta_p_SR,
+ FW2Data& fW_2) const;
+
+ void dEval(BiotData const& biot_data,
+ CapillaryPressureData const pCap,
+ ConstituentDensityData const& constituent_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PorosityDerivativeData const& porosity_d_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ SolidCompressibilityData const& beta_p_SR,
+ FW2DerivativeData& dfW_2) const;
+};
+
+struct FW3aData
+{
+ double a = nan;
+};
+
+struct FW3aDerivativeData
+{
+ double dp_GR = nan;
+ double dp_cap = nan;
+ double dT = nan;
+};
+
+struct FW3aModel
+{
+ void eval(
+ double const dt,
+ ConstituentDensityData const& constituent_density_data,
+ PrevState<ConstituentDensityData> const& constituent_density_data_prev,
+ PrevState<PureLiquidDensityData> const& rho_W_LR_prev,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ FW3aData& fW_3a) const;
+
+ void dEval(
+ double const dt,
+ ConstituentDensityData const& constituent_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PrevState<ConstituentDensityData> const& constituent_density_data_prev,
+ PrevState<PureLiquidDensityData> const& rho_W_LR_prev,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ FW3aDerivativeData& dfW_3a) const;
+};
+
+template <int DisplacementDim>
+struct FW4LWpGData
+{
+ GlobalDimMatrix<DisplacementDim> L;
+};
+
+template <int DisplacementDim>
+struct FW4LWpGDerivativeData
+{
+ GlobalDimMatrix<DisplacementDim> dp_GR;
+ GlobalDimMatrix<DisplacementDim> dp_cap;
+ GlobalDimMatrix<DisplacementDim> dT;
+};
+
+template <int DisplacementDim>
+struct FW4LWpGModel
+{
+ void eval(AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FW4LWpGData<DisplacementDim>& fW_4_LWpG) const;
+
+ void dEval(ConstituentDensityData const& constituent_density_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ PhaseTransitionData const& phase_transition_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ ViscosityData const& viscosity_data,
+ FW4LWpGDerivativeData<DisplacementDim>& dfW_4_LWpG) const;
+};
+
+extern template struct FW4LWpGModel<2>;
+extern template struct FW4LWpGModel<3>;
+
+template <int DisplacementDim>
+struct FW4LWpCData
+{
+ GlobalDimMatrix<DisplacementDim> L;
+};
+
+template <int DisplacementDim>
+struct FW4LWpCDerivativeData
+{
+ GlobalDimMatrix<DisplacementDim> dp_GR;
+ GlobalDimMatrix<DisplacementDim> dp_cap;
+ GlobalDimMatrix<DisplacementDim> dT;
+};
+
+template <int DisplacementDim>
+struct FW4LWpCModel
+{
+ void eval(AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FW4LWpCData<DisplacementDim>& fW_4_LWpC) const;
+
+ void dEval(AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ ViscosityData const& viscosity_data,
+ FW4LWpCDerivativeData<DisplacementDim>& dfW_4_LWpC) const;
+};
+
+extern template struct FW4LWpCModel<2>;
+extern template struct FW4LWpCModel<3>;
+
+template <int DisplacementDim>
+struct FW4LWTData
+{
+ GlobalDimMatrix<DisplacementDim> L;
+};
+
+template <int DisplacementDim>
+struct FW4LWTModel
+{
+ void eval(FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FW4LWTData<DisplacementDim>& fW_4_LWT) const;
+};
+
+struct FW4MWpGData
+{
+ double m = nan;
+};
+
+struct FW4MWpGModel
+{
+ void eval(BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const& beta_p_SR,
+ FW4MWpGData& fW_4_MWpG) const;
+};
+
+struct FW4MWpCData
+{
+ double m = nan;
+ double ml = nan;
+};
+
+struct FW4MWpCModel
+{
+ void eval(BiotData const& biot_data,
+ CapillaryPressureData const pCap,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ PrevState<SaturationData> const& S_L_data_prev,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SolidCompressibilityData const& beta_p_SR,
+ FW4MWpCData& fW_4_MWpC) const;
+};
+
+struct FW4MWTData
+{
+ double m = nan;
+};
+
+template <int DisplacementDim>
+struct FW4MWTModel
+{
+ void eval(
+ BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ FW4MWTData& fW_4_MWT) const;
+};
+
+extern template struct FW4MWTModel<2>;
+extern template struct FW4MWTModel<3>;
+
+struct FW4MWuData
+{
+ double m = nan;
+};
+
+struct FW4MWuModel
+{
+ void eval(BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ FW4MWuData& fW_4_MWu) const;
+};
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 2f14daaaeebb1d46f078dd6b93b23d245cc5ab5b..154248a716b36a6b3576d6aab655ab79d156e17d 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -152,16 +152,19 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
double const T = NT.dot(temperature);
double const T_prev = NT.dot(temperature_prev);
- double const pGR = Np.dot(gas_pressure);
- double const pCap = Np.dot(capillary_pressure);
ConstitutiveRelations::TemperatureData const T_data{T, T_prev};
- ConstitutiveRelations::GasPressureData const pGR_data{pGR};
- ConstitutiveRelations::CapillaryPressureData const pCap_data{pCap};
+ ConstitutiveRelations::GasPressureData const pGR_data{
+ Np.dot(gas_pressure)};
+ ConstitutiveRelations::CapillaryPressureData const pCap_data{
+ Np.dot(capillary_pressure)};
ConstitutiveRelations::ReferenceTemperatureData const T0{
this->process_data_.reference_temperature(t, pos)[0]};
- GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
- GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
- GlobalDimVectorType const gradT = gradNp * temperature;
+ ConstitutiveRelations::GasPressureGradientData<DisplacementDim> const
+ grad_p_GR{gradNp * gas_pressure};
+ ConstitutiveRelations::CapillaryPressureGradientData<
+ DisplacementDim> const grad_p_cap{gradNp * capillary_pressure};
+ ConstitutiveRelations::TemperatureGradientData<DisplacementDim> const
+ grad_T{gradNp * temperature};
// medium properties
models.elastic_tangent_stiffness_model.eval({pos, t, dt}, T_data,
@@ -175,10 +178,9 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
typename BMatricesType::BMatrixType>(
gradNu, Nu, x_coord, this->is_axially_symmetric_);
- auto& eps = ip_out.eps_data.eps;
- eps.noalias() = Bu * displacement;
+ ip_out.eps_data.eps.noalias() = Bu * displacement;
models.S_L_model.eval({pos, t, dt}, media_data, pCap_data,
- current_state.S_L_data, ip_cv.dS_L_dp_cap);
+ current_state.S_L_data);
models.chi_S_L_model.eval({pos, t, dt}, media_data,
current_state.S_L_data, ip_cv.chi_S_L);
@@ -223,7 +225,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
models.phase_transition_model.eval(
{pos, t, dt}, media_data, pGR_data, pCap_data, T_data,
- current_state.rho_W_LR, ip_out.enthalpy_data,
+ current_state.rho_W_LR, ip_out.fluid_enthalpy_data,
ip_out.mass_mole_fractions_data, ip_out.fluid_density_data,
ip_out.vapour_pressure_data, current_state.constituent_density_data,
ip_cv.phase_transition_data);
@@ -237,488 +239,496 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
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);
+ ip_out.porosity_data);
- models.solid_density_model.eval(
- {pos, t, dt}, media_data, T_data,
+ models.solid_density_model.eval({pos, t, dt}, media_data, T_data,
#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- ip_cv.biot_data, ip_out.eps_data, ip_cv.s_therm_exp_data,
+ ip_cv.biot_data, ip_out.eps_data,
+ ip_cv.s_therm_exp_data,
#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- ip_out.solid_density_data, ip_cv.solid_density_d_data);
+ ip_out.solid_density_data);
models.solid_heat_capacity_model.eval({pos, t, dt}, media_data, T_data,
ip_cv.solid_heat_capacity_data);
models.thermal_conductivity_model.eval(
{pos, t, dt}, media_data, T_data, ip_out.porosity_data,
- ip_cv.porosity_d_data, current_state.S_L_data, ip_cv.dS_L_dp_cap,
- ip_cv.thermal_conductivity_data);
-
- auto const& c = ip_cv.phase_transition_data;
-
- auto const phi_L =
- 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;
-
- ip_out.enthalpy_data.h_S = ip_cv.solid_heat_capacity_data() * T;
- auto const u_S = ip_out.enthalpy_data.h_S;
-
- current_state.internal_energy_data() =
- phi_G * ip_out.fluid_density_data.rho_GR * c.uG +
- phi_L * ip_out.fluid_density_data.rho_LR * c.uL +
- phi_S * ip_out.solid_density_data.rho_SR * u_S;
-
- ip_cv.effective_volumetric_enthalpy_data.rho_h_eff =
- phi_G * ip_out.fluid_density_data.rho_GR *
- ip_out.enthalpy_data.h_G +
- phi_L * ip_out.fluid_density_data.rho_LR *
- ip_out.enthalpy_data.h_L +
- phi_S * ip_out.solid_density_data.rho_SR * ip_out.enthalpy_data.h_S;
-
- // for variable output
- auto const xmCL = 1. - ip_out.mass_mole_fractions_data.xmWL;
-
- const GlobalDimVectorType gradxmWG = c.dxmWG_dpGR * gradpGR +
- c.dxmWG_dpCap * gradpCap +
- c.dxmWG_dT * gradT;
- const GlobalDimVectorType gradxmCG = -gradxmWG;
-
- const GlobalDimVectorType gradxmWL = c.dxmWL_dpGR * gradpGR +
- c.dxmWL_dpCap * gradpCap +
- c.dxmWL_dT * gradT;
- const GlobalDimVectorType gradxmCL = -gradxmWL;
-
- // Todo: factor -phiAlpha / xmZetaAlpha * DZetaAlpha can be evaluated in
- // the respective phase transition model, here only the multiplication
- // with the gradient of the mass fractions should take place.
-
- ip_out.diffusion_velocity_data.d_CG =
- ip_out.mass_mole_fractions_data.xmCG == 0.
- ? 0. * gradxmCG // Keep d_CG's dimension and prevent
- // division by zero
- : -phi_G / ip_out.mass_mole_fractions_data.xmCG *
- c.diffusion_coefficient_vapour * gradxmCG;
-
- ip_out.diffusion_velocity_data.d_WG =
- ip_out.mass_mole_fractions_data.xmCG == 1.
- ? 0. * gradxmWG // Keep d_WG's dimension and prevent
- // division by zero
- : -phi_G / (1 - ip_out.mass_mole_fractions_data.xmCG) *
- c.diffusion_coefficient_vapour * gradxmWG;
-
- ip_out.diffusion_velocity_data.d_CL =
- xmCL == 0.
- ? 0. * gradxmCL // Keep d_CL's dimension and
- // prevent division by zero
- : -phi_L / xmCL * c.diffusion_coefficient_solute * gradxmCL;
-
- ip_out.diffusion_velocity_data.d_WL =
- ip_out.mass_mole_fractions_data.xmWL == 0.
- ? 0. * gradxmWL // Keep d_WG's dimension and prevent
- // division by zero
- : -phi_L / ip_out.mass_mole_fractions_data.xmWL *
- c.diffusion_coefficient_solute * gradxmWL;
-
- // ---------------------------------------------------------------------
- // Derivatives for Jacobian
- // ---------------------------------------------------------------------
- ip_cv.effective_volumetric_internal_energy_d_data.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 +
- phi_L * c.drho_LR_dT * c.uL +
- phi_L * ip_out.fluid_density_data.rho_LR * c.du_L_dT +
- phi_S * ip_cv.solid_density_d_data.drho_SR_dT * u_S +
- phi_S * ip_out.solid_density_data.rho_SR *
- ip_cv.solid_heat_capacity_data() -
- ip_cv.porosity_d_data.dphi_dT * ip_out.solid_density_data.rho_SR *
- u_S;
-
- // ds_L_dp_GR = 0;
- // ds_G_dp_GR = -ds_L_dp_GR;
- double const ds_G_dp_cap = -ip_cv.dS_L_dp_cap();
-
- // dphi_G_dp_GR = -ds_L_dp_GR * ip_out.porosity_data.phi = 0;
- double const dphi_G_dp_cap =
- -ip_cv.dS_L_dp_cap() * ip_out.porosity_data.phi;
- // dphi_L_dp_GR = ds_L_dp_GR * ip_out.porosity_data.phi = 0;
- double const dphi_L_dp_cap =
- ip_cv.dS_L_dp_cap() * ip_out.porosity_data.phi;
-
- // From p_LR = p_GR - p_cap it follows for
- // drho_LR/dp_GR = drho_LR/dp_LR * dp_LR/dp_GR
- // = drho_LR/dp_LR * (dp_GR/dp_GR - dp_cap/dp_GR)
- // = drho_LR/dp_LR * (1 - 0)
- double const drho_LR_dp_GR = c.drho_LR_dp_LR;
- double const drho_LR_dp_cap = -c.drho_LR_dp_LR;
- // drho_GR_dp_cap = 0;
-
- ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_GR =
- /*(dphi_G_dp_GR = 0) * ip_out.fluid_density_data.rho_GR *
- ip_out.enthalpy_data.h_G +*/
- phi_G * c.drho_GR_dp_GR * ip_out.enthalpy_data.h_G +
- /*(dphi_L_dp_GR = 0) * ip_out.fluid_density_data.rho_LR *
- ip_out.enthalpy_data.h_L +*/
- phi_L * drho_LR_dp_GR * ip_out.enthalpy_data.h_L;
- ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_cap =
- dphi_G_dp_cap * ip_out.fluid_density_data.rho_GR *
- ip_out.enthalpy_data.h_G +
- /*phi_G * (drho_GR_dp_cap = 0) * ip_out.enthalpy_data.h_G +*/
- dphi_L_dp_cap * ip_out.fluid_density_data.rho_LR *
- ip_out.enthalpy_data.h_L +
- phi_L * drho_LR_dp_cap * ip_out.enthalpy_data.h_L;
-
- // TODO (naumov) Extend for temperature dependent porosities.
- constexpr double dphi_G_dT = 0;
- constexpr double dphi_L_dT = 0;
- ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dT =
- dphi_G_dT * ip_out.fluid_density_data.rho_GR *
- ip_out.enthalpy_data.h_G +
- phi_G * c.drho_GR_dT * ip_out.enthalpy_data.h_G +
- phi_G * ip_out.fluid_density_data.rho_GR * c.dh_G_dT +
- dphi_L_dT * ip_out.fluid_density_data.rho_LR *
- ip_out.enthalpy_data.h_L +
- phi_L * c.drho_LR_dT * ip_out.enthalpy_data.h_L +
- phi_L * ip_out.fluid_density_data.rho_LR * c.dh_L_dT -
- ip_cv.porosity_d_data.dphi_dT * ip_out.solid_density_data.rho_SR *
- ip_out.enthalpy_data.h_S +
- phi_S * ip_cv.solid_density_d_data.drho_SR_dT *
- ip_out.enthalpy_data.h_S +
- phi_S * ip_out.solid_density_data.rho_SR *
- ip_cv.solid_heat_capacity_data();
-
- ip_cv.effective_volumetric_internal_energy_d_data.drho_u_eff_dp_GR =
- /*(dphi_G_dp_GR = 0) * ip_out.fluid_density_data.rho_GR * c.uG +*/
- phi_G * c.drho_GR_dp_GR * c.uG +
- phi_G * ip_out.fluid_density_data.rho_GR * c.du_G_dp_GR +
- /*(dphi_L_dp_GR = 0) * ip_out.fluid_density_data.rho_LR * c.uL +*/
- phi_L * drho_LR_dp_GR * c.uL +
- phi_L * ip_out.fluid_density_data.rho_LR * c.du_L_dp_GR;
-
- ip_cv.effective_volumetric_internal_energy_d_data.drho_u_eff_dp_cap =
- dphi_G_dp_cap * ip_out.fluid_density_data.rho_GR * c.uG +
- /*phi_G * (drho_GR_dp_cap = 0) * c.uG +*/
- dphi_L_dp_cap * ip_out.fluid_density_data.rho_LR * c.uL +
- phi_L * drho_LR_dp_cap * c.uL +
- phi_L * ip_out.fluid_density_data.rho_LR * c.du_L_dp_cap;
+ current_state.S_L_data, ip_cv.thermal_conductivity_data);
+
+ models.advection_model.eval(current_state.constituent_density_data,
+ ip_out.permeability_data,
+ current_state.rho_W_LR,
+ ip_cv.viscosity_data,
+ ip_cv.advection_data);
+
+ models.gravity_model.eval(
+ ip_out.fluid_density_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_out.solid_density_data,
+ ConstitutiveRelations::SpecificBodyForceData<DisplacementDim>{
+ this->process_data_.specific_body_force},
+ ip_cv.volumetric_body_force);
+
+ models.diffusion_velocity_model.eval(grad_p_cap,
+ grad_p_GR,
+ ip_out.mass_mole_fractions_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ grad_T,
+ ip_out.diffusion_velocity_data);
+
+ models.solid_enthalpy_model.eval(ip_cv.solid_heat_capacity_data, T_data,
+ ip_out.solid_enthalpy_data);
+
+ models.internal_energy_model.eval(ip_out.fluid_density_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_out.solid_density_data,
+ ip_out.solid_enthalpy_data,
+ current_state.internal_energy_data);
+
+ models.effective_volumetric_enthalpy_model.eval(
+ ip_out.fluid_density_data,
+ ip_out.fluid_enthalpy_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_out.solid_density_data,
+ ip_out.solid_enthalpy_data,
+ ip_cv.effective_volumetric_enthalpy_data);
+
+ models.fC_1_model.eval(ip_cv.advection_data, ip_out.fluid_density_data,
+ ip_cv.fC_1);
- auto const& b = this->process_data_.specific_body_force;
- GlobalDimMatrixType const k_over_mu_G =
- ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
- ip_cv.viscosity_data.mu_GR;
- GlobalDimMatrixType const k_over_mu_L =
- ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
- ip_cv.viscosity_data.mu_LR;
-
- // dk_over_mu_G_dp_GR = ip_out.permeability_data.Ki *
- // ip_out.permeability_data.dk_rel_G_dS_L *
- // (ds_L_dp_GR = 0) /
- // ip_cv.viscosity_data.mu_GR = 0;
- // dk_over_mu_L_dp_GR = ip_out.permeability_data.Ki *
- // ip_out.permeability_data.dk_rel_L_dS_L *
- // (ds_L_dp_GR = 0) /
- // ip_cv.viscosity_data.mu_LR = 0;
- ip_cv.dk_over_mu_G_dp_cap = ip_out.permeability_data.Ki *
- ip_out.permeability_data.dk_rel_G_dS_L *
- ip_cv.dS_L_dp_cap() /
- ip_cv.viscosity_data.mu_GR;
- ip_cv.dk_over_mu_L_dp_cap = ip_out.permeability_data.Ki *
- ip_out.permeability_data.dk_rel_L_dS_L *
- ip_cv.dS_L_dp_cap() /
- ip_cv.viscosity_data.mu_LR;
-
- ip_out.darcy_velocity_data.w_GS =
- k_over_mu_G * ip_out.fluid_density_data.rho_GR * b -
- k_over_mu_G * gradpGR;
- ip_out.darcy_velocity_data.w_LS =
- k_over_mu_L * gradpCap +
- k_over_mu_L * ip_out.fluid_density_data.rho_LR * b -
- k_over_mu_L * gradpGR;
-
- ip_cv.drho_GR_h_w_eff_dp_GR_Npart =
- c.drho_GR_dp_GR * ip_out.enthalpy_data.h_G *
- ip_out.darcy_velocity_data.w_GS +
- ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G *
- k_over_mu_G * c.drho_GR_dp_GR * b;
- ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart =
- ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G *
- k_over_mu_G -
- ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L *
- k_over_mu_L;
-
- ip_cv.drho_LR_h_w_eff_dp_cap_Npart =
- -drho_LR_dp_cap * ip_out.enthalpy_data.h_L *
- ip_out.darcy_velocity_data.w_LS -
- ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L *
- k_over_mu_L * drho_LR_dp_cap * b;
- ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart =
- // TODO (naumov) why the minus sign??????
- ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L *
- k_over_mu_L;
-
- ip_cv.drho_GR_h_w_eff_dT =
- c.drho_GR_dT * ip_out.enthalpy_data.h_G *
- ip_out.darcy_velocity_data.w_GS +
- ip_out.fluid_density_data.rho_GR * c.dh_G_dT *
- ip_out.darcy_velocity_data.w_GS +
- c.drho_LR_dT * ip_out.enthalpy_data.h_L *
- ip_out.darcy_velocity_data.w_LS +
- ip_out.fluid_density_data.rho_LR * c.dh_L_dT *
- ip_out.darcy_velocity_data.w_LS;
- // TODO (naumov) + k_over_mu_G * drho_GR_dT * b + k_over_mu_L *
- // drho_LR_dT * b
-
- // Derivatives of s_G * rho_C_GR_dot + s_L * rho_C_LR_dot abbreviated
- // here with S_rho_C_eff.
- double const s_L = current_state.S_L_data.S_L;
- double const s_G = 1. - s_L;
- double const rho_C_FR =
- s_G * current_state.constituent_density_data.rho_C_GR +
- s_L * current_state.constituent_density_data.rho_C_LR;
- double const rho_W_FR =
- s_G * current_state.constituent_density_data.rho_W_GR +
- s_L * current_state.rho_W_LR();
- // TODO (naumov) Extend for partially saturated media.
- constexpr double drho_C_GR_dp_cap = 0;
- if (dt == 0.)
+ if (!this->process_data_.apply_mass_lumping)
{
- ip_cv.dfC_3a_dp_GR = 0.;
- ip_cv.dfC_3a_dp_cap = 0.;
- ip_cv.dfC_3a_dT = 0.;
+ models.fC_2a_model.eval(ip_cv.biot_data,
+ pCap_data,
+ current_state.constituent_density_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_cv.beta_p_SR,
+ ip_cv.fC_2a);
}
- else
+ models.fC_3a_model.eval(dt,
+ current_state.constituent_density_data,
+ prev_state.constituent_density_data,
+ current_state.S_L_data,
+ ip_cv.fC_3a);
+
+ models.fC_4_LCpG_model.eval(ip_cv.advection_data,
+ ip_out.fluid_density_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_cv.fC_4_LCpG);
+
+ models.fC_4_LCpC_model.eval(ip_cv.advection_data,
+ ip_out.fluid_density_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_cv.fC_4_LCpC);
+
+ models.fC_4_LCT_model.eval(ip_out.fluid_density_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_cv.fC_4_LCT);
+
+ models.fC_4_MCpG_model.eval(ip_cv.biot_data,
+ current_state.constituent_density_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_cv.beta_p_SR,
+ ip_cv.fC_4_MCpG);
+
+ models.fC_4_MCpC_model.eval(ip_cv.biot_data,
+ pCap_data,
+ current_state.constituent_density_data,
+ ip_out.porosity_data,
+ prev_state.S_L_data,
+ current_state.S_L_data,
+ ip_cv.beta_p_SR,
+ ip_cv.fC_4_MCpC);
+
+ models.fC_4_MCT_model.eval(ip_cv.biot_data,
+ current_state.constituent_density_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_cv.s_therm_exp_data,
+ ip_cv.fC_4_MCT);
+
+ models.fC_4_MCu_model.eval(ip_cv.biot_data,
+ current_state.constituent_density_data,
+ current_state.S_L_data,
+ ip_cv.fC_4_MCu);
+
+ models.fW_1_model.eval(ip_cv.advection_data, ip_out.fluid_density_data,
+ ip_cv.fW_1);
+
+ if (!this->process_data_.apply_mass_lumping)
{
- double const rho_C_GR_dot =
- (current_state.constituent_density_data.rho_C_GR -
- prev_state.constituent_density_data->rho_C_GR) /
- dt;
- double const rho_C_LR_dot =
- (current_state.constituent_density_data.rho_C_LR -
- prev_state.constituent_density_data->rho_C_LR) /
- dt;
- ip_cv.dfC_3a_dp_GR =
- /*(ds_G_dp_GR = 0) * rho_C_GR_dot +*/ s_G * c.drho_C_GR_dp_GR /
- dt +
- /*(ds_L_dp_GR = 0) * rho_C_LR_dot +*/ s_L * c.drho_C_LR_dp_GR /
- dt;
- ip_cv.dfC_3a_dp_cap = ds_G_dp_cap * rho_C_GR_dot +
- s_G * drho_C_GR_dp_cap / dt +
- ip_cv.dS_L_dp_cap() * rho_C_LR_dot -
- s_L * c.drho_C_LR_dp_LR / dt;
- ip_cv.dfC_3a_dT =
- s_G * c.drho_C_GR_dT / dt + s_L * c.drho_C_LR_dT / dt;
+ models.fW_2_model.eval(ip_cv.biot_data,
+ pCap_data,
+ current_state.constituent_density_data,
+ ip_out.porosity_data,
+ current_state.rho_W_LR,
+ current_state.S_L_data,
+ ip_cv.beta_p_SR,
+ ip_cv.fW_2);
}
+ models.fW_3a_model.eval(dt,
+ current_state.constituent_density_data,
+ prev_state.constituent_density_data,
+ prev_state.rho_W_LR,
+ current_state.rho_W_LR,
+ current_state.S_L_data,
+ ip_cv.fW_3a);
+
+ models.fW_4_LWpG_model.eval(ip_cv.advection_data,
+ ip_out.fluid_density_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_cv.fW_4_LWpG);
+
+ models.fW_4_LWpC_model.eval(ip_cv.advection_data,
+ ip_out.fluid_density_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_cv.fW_4_LWpC);
+
+ models.fW_4_LWT_model.eval(ip_out.fluid_density_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_cv.fW_4_LWT);
+
+ models.fW_4_MWpG_model.eval(ip_cv.biot_data,
+ current_state.constituent_density_data,
+ ip_out.porosity_data,
+ current_state.rho_W_LR,
+ current_state.S_L_data,
+ ip_cv.beta_p_SR,
+ ip_cv.fW_4_MWpG);
+
+ models.fW_4_MWpC_model.eval(ip_cv.biot_data,
+ pCap_data,
+ current_state.constituent_density_data,
+ ip_out.porosity_data,
+ prev_state.S_L_data,
+ current_state.rho_W_LR,
+ current_state.S_L_data,
+ ip_cv.beta_p_SR,
+ ip_cv.fW_4_MWpC);
+
+ models.fW_4_MWT_model.eval(ip_cv.biot_data,
+ current_state.constituent_density_data,
+ ip_out.porosity_data,
+ current_state.rho_W_LR,
+ current_state.S_L_data,
+ ip_cv.s_therm_exp_data,
+ ip_cv.fW_4_MWT);
+
+ models.fW_4_MWu_model.eval(ip_cv.biot_data,
+ current_state.constituent_density_data,
+ current_state.rho_W_LR,
+ current_state.S_L_data,
+ ip_cv.fW_4_MWu);
+
+ models.fT_1_model.eval(dt,
+ current_state.internal_energy_data,
+ prev_state.internal_energy_data,
+ ip_cv.fT_1);
+
+ // ---------------------------------------------------------------------
+ // Derivatives for Jacobian
+ // ---------------------------------------------------------------------
+
+ models.darcy_velocity_model.eval(
+ grad_p_cap,
+ ip_out.fluid_density_data,
+ grad_p_GR,
+ ip_out.permeability_data,
+ ConstitutiveRelations::SpecificBodyForceData<DisplacementDim>{
+ this->process_data_.specific_body_force},
+ ip_cv.viscosity_data,
+ ip_out.darcy_velocity_data);
+
+ models.fT_2_model.eval(ip_out.darcy_velocity_data,
+ ip_out.fluid_density_data,
+ ip_out.fluid_enthalpy_data,
+ ip_cv.fT_2);
+
+ models.fT_3_model.eval(
+ current_state.constituent_density_data,
+ ip_out.darcy_velocity_data,
+ ip_out.diffusion_velocity_data,
+ ip_out.fluid_density_data,
+ ip_cv.phase_transition_data,
+ ConstitutiveRelations::SpecificBodyForceData<DisplacementDim>{
+ this->process_data_.specific_body_force},
+ ip_cv.fT_3);
+
+ models.fu_2_KupC_model.eval(ip_cv.biot_data, ip_cv.chi_S_L,
+ ip_cv.fu_2_KupC);
+ }
- double const drho_C_FR_dp_GR =
- /*(ds_G_dp_GR = 0) * current_state.constituent_density_data.rho_C_GR
- +*/
- s_G * c.drho_C_GR_dp_GR +
- /*(ds_L_dp_GR = 0) * current_state.constituent_density_data.rho_C_LR
- +*/
- s_L * c.drho_C_LR_dp_GR;
- ip_cv.dfC_4_MCpG_dp_GR =
- drho_C_FR_dp_GR * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
- ip_cv.beta_p_SR();
-
- double const drho_C_FR_dT = s_G * c.drho_C_GR_dT + s_L * c.drho_C_LR_dT;
- ip_cv.dfC_4_MCpG_dT =
- drho_C_FR_dT * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
- ip_cv.beta_p_SR() -
- rho_C_FR * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR();
-
- ip_cv.dfC_4_MCT_dT =
- drho_C_FR_dT * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
- ip_cv.s_therm_exp_data.beta_T_SR
+ return {ip_constitutive_data, ip_constitutive_variables};
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ int DisplacementDim>
+std::vector<ConstitutiveRelations::DerivativesData<DisplacementDim>>
+TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
+ DisplacementDim>::
+ updateConstitutiveVariablesDerivatives(
+ Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev,
+ double const t, double const dt,
+ std::vector<
+ ConstitutiveRelations::ConstitutiveData<DisplacementDim>> const&
+ ip_constitutive_data,
+ std::vector<
+ ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>> const&
+ ip_constitutive_variables,
+ ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const&
+ models)
+{
+ [[maybe_unused]] auto const matrix_size =
+ gas_pressure_size + capillary_pressure_size + temperature_size +
+ displacement_size;
+
+ assert(local_x.size() == matrix_size);
+
+ auto const temperature =
+ local_x.template segment<temperature_size>(temperature_index);
+ auto const temperature_prev =
+ local_x_prev.template segment<temperature_size>(temperature_index);
+
+ auto const capillary_pressure =
+ local_x.template segment<capillary_pressure_size>(
+ capillary_pressure_index);
+
+ auto const& medium =
+ *this->process_data_.media_map.getMedium(this->element_.getID());
+ ConstitutiveRelations::MediaData media_data{medium};
+
+ unsigned const n_integration_points =
+ this->integration_method_.getNumberOfPoints();
+
+ std::vector<ConstitutiveRelations::DerivativesData<DisplacementDim>>
+ ip_d_data(n_integration_points);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& ip_data = _ip_data[ip];
+ auto& ip_dd = ip_d_data[ip];
+ auto const& ip_cd = ip_constitutive_data[ip];
+ auto const& ip_cv = ip_constitutive_variables[ip];
+ auto const& ip_out = this->output_data_[ip];
+ auto const& current_state = this->current_states_[ip];
+ auto const& prev_state = this->prev_states_[ip];
+
+ auto const& Nu = ip_data.N_u;
+ auto const& Np = ip_data.N_p;
+ auto const& NT = Np;
+
+ ParameterLib::SpatialPosition const pos{
+ std::nullopt, this->element_.getID(), ip,
+ MathLib::Point3d(
+ NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement>(
+ this->element_, Nu))};
+
+ double const T = NT.dot(temperature);
+ double const T_prev = NT.dot(temperature_prev);
+ ConstitutiveRelations::TemperatureData const T_data{T, T_prev};
+ ConstitutiveRelations::CapillaryPressureData const pCap_data{
+ Np.dot(capillary_pressure)};
+
+ models.S_L_model.dEval({pos, t, dt}, media_data, pCap_data,
+ ip_dd.dS_L_dp_cap);
+
+ models.advection_model.dEval(current_state.constituent_density_data,
+ ip_out.permeability_data,
+ ip_cv.viscosity_data,
+ ip_dd.dS_L_dp_cap,
+ ip_cv.phase_transition_data,
+ ip_dd.advection_d_data);
+
+ models.porosity_model.dEval(
+ {pos, t, dt}, media_data, ip_out.porosity_data, ip_dd.dS_L_dp_cap,
#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- + rho_C_FR * (ip_cv.biot_data() - ip_cv.porosity_d_data.dphi_dT) *
- ip_cv.s_therm_exp_data.beta_T_SR
-#endif
- ;
-
- ip_cv.dfC_4_MCu_dT = drho_C_FR_dT * ip_cv.biot_data();
-
- ip_cv.dfC_2a_dp_GR =
- -ip_out.porosity_data.phi * c.drho_C_GR_dp_GR -
- drho_C_FR_dp_GR * pCap *
- (ip_cv.biot_data() - ip_out.porosity_data.phi) *
- ip_cv.beta_p_SR();
-
- double const drho_C_FR_dp_cap =
- ds_G_dp_cap * current_state.constituent_density_data.rho_C_GR +
- s_G * drho_C_GR_dp_cap +
- ip_cv.dS_L_dp_cap() *
- current_state.constituent_density_data.rho_C_LR -
- s_L * c.drho_C_LR_dp_LR;
-
- ip_cv.dfC_2a_dp_cap =
- ip_out.porosity_data.phi * (-c.drho_C_LR_dp_LR - drho_C_GR_dp_cap) -
- drho_C_FR_dp_cap * pCap *
- (ip_cv.biot_data() - ip_out.porosity_data.phi) *
- ip_cv.beta_p_SR() +
- rho_C_FR * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
- ip_cv.beta_p_SR();
-
- ip_cv.dfC_2a_dT =
- ip_cv.porosity_d_data.dphi_dT *
- (current_state.constituent_density_data.rho_C_LR -
- current_state.constituent_density_data.rho_C_GR) +
- ip_out.porosity_data.phi * (c.drho_C_LR_dT - c.drho_C_GR_dT) -
- drho_C_FR_dT * pCap *
- (ip_cv.biot_data() - ip_out.porosity_data.phi) *
- ip_cv.beta_p_SR() +
- rho_C_FR * pCap * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR();
-
- ip_cv.dadvection_C_dp_GR = c.drho_C_GR_dp_GR * k_over_mu_G
- // + rhoCGR * (dk_over_mu_G_dp_GR = 0)
- // + rhoCLR * (dk_over_mu_L_dp_GR = 0)
- + c.drho_C_LR_dp_GR * k_over_mu_L;
-
- ip_cv.dadvection_C_dp_cap =
- //(drho_C_GR_dp_cap = 0) * k_over_mu_G
- current_state.constituent_density_data.rho_C_GR *
- ip_cv.dk_over_mu_G_dp_cap +
- (-c.drho_C_LR_dp_LR) * k_over_mu_L +
- current_state.constituent_density_data.rho_C_LR *
- ip_cv.dk_over_mu_L_dp_cap;
-
- ip_cv.dfC_4_LCpG_dT =
- c.drho_C_GR_dT * k_over_mu_G + c.drho_C_LR_dT * k_over_mu_L
- // + ip_cv.ddiffusion_C_p_dT TODO (naumov)
- ;
-
- double const drho_W_FR_dp_GR =
- /*(ds_G_dp_GR = 0) * current_state.constituent_density_data.rho_W_GR
- +*/
- s_G * c.drho_W_GR_dp_GR +
- /*(ds_L_dp_GR = 0) * current_state.rho_W_LR() +*/
- s_L * c.drho_W_LR_dp_GR;
- double const drho_W_FR_dp_cap =
- ds_G_dp_cap * current_state.constituent_density_data.rho_W_GR +
- s_G * c.drho_W_GR_dp_cap +
- ip_cv.dS_L_dp_cap() * current_state.rho_W_LR() -
- s_L * c.drho_W_LR_dp_LR;
- double const drho_W_FR_dT = s_G * c.drho_W_GR_dT + s_L * c.drho_W_LR_dT;
-
- ip_cv.dfW_2a_dp_GR =
- ip_out.porosity_data.phi * (c.drho_W_LR_dp_GR - c.drho_W_GR_dp_GR);
- ip_cv.dfW_2b_dp_GR = drho_W_FR_dp_GR * pCap *
- (ip_cv.biot_data() - ip_out.porosity_data.phi) *
- ip_cv.beta_p_SR();
- ip_cv.dfW_2a_dp_cap = ip_out.porosity_data.phi *
- (-c.drho_W_LR_dp_LR - c.drho_W_GR_dp_cap);
- ip_cv.dfW_2b_dp_cap =
- drho_W_FR_dp_cap * pCap *
- (ip_cv.biot_data() - ip_out.porosity_data.phi) *
- ip_cv.beta_p_SR() +
- rho_W_FR * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
- ip_cv.beta_p_SR();
-
- ip_cv.dfW_2a_dT =
- ip_cv.porosity_d_data.dphi_dT *
- (current_state.rho_W_LR() -
- current_state.constituent_density_data.rho_W_GR) +
- ip_out.porosity_data.phi * (c.drho_W_LR_dT - c.drho_W_GR_dT);
- ip_cv.dfW_2b_dT =
- drho_W_FR_dT * pCap *
- (ip_cv.biot_data() - ip_out.porosity_data.phi) *
- ip_cv.beta_p_SR() -
- rho_W_FR * pCap * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR();
-
- if (dt == 0.)
+ ip_cv.biot_data, ip_out.eps_data, ip_cv.s_therm_exp_data,
+#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
+ ip_dd.porosity_d_data);
+
+ models.thermal_conductivity_model.dEval(
+ {pos, t, dt}, media_data, T_data, ip_out.porosity_data,
+ ip_dd.porosity_d_data, current_state.S_L_data,
+ ip_dd.thermal_conductivity_d_data);
+
+ models.solid_density_model.dEval({pos, t, dt}, media_data, T_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_dd.solid_density_d_data);
+
+ models.internal_energy_model.dEval(
+ ip_out.fluid_density_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ ip_dd.porosity_d_data,
+ current_state.S_L_data,
+ ip_out.solid_density_data,
+ ip_dd.solid_density_d_data,
+ ip_out.solid_enthalpy_data,
+ ip_cv.solid_heat_capacity_data,
+ ip_dd.effective_volumetric_internal_energy_d_data);
+
+ models.effective_volumetric_enthalpy_model.dEval(
+ ip_out.fluid_density_data,
+ ip_out.fluid_enthalpy_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ ip_dd.porosity_d_data,
+ current_state.S_L_data,
+ ip_out.solid_density_data,
+ ip_dd.solid_density_d_data,
+ ip_out.solid_enthalpy_data,
+ ip_cv.solid_heat_capacity_data,
+ ip_dd.effective_volumetric_enthalpy_d_data);
+ if (!this->process_data_.apply_mass_lumping)
{
- ip_cv.dfW_3a_dp_GR = 0.;
- ip_cv.dfW_3a_dp_cap = 0.;
- ip_cv.dfW_3a_dT = 0.;
+ models.fC_2a_model.dEval(ip_cv.biot_data,
+ pCap_data,
+ current_state.constituent_density_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ ip_dd.porosity_d_data,
+ current_state.S_L_data,
+ ip_dd.dS_L_dp_cap,
+ ip_cv.beta_p_SR,
+ ip_dd.dfC_2a);
}
- else
+ models.fC_3a_model.dEval(dt,
+ current_state.constituent_density_data,
+ prev_state.constituent_density_data,
+ ip_cv.phase_transition_data,
+ current_state.S_L_data,
+ ip_dd.dS_L_dp_cap,
+ ip_dd.dfC_3a);
+
+ models.fC_4_LCpG_model.dEval(ip_out.permeability_data,
+ ip_cv.viscosity_data,
+ ip_cv.phase_transition_data,
+ ip_dd.advection_d_data,
+ ip_dd.dfC_4_LCpG);
+
+ models.fC_4_LCpC_model.dEval(current_state.constituent_density_data,
+ ip_out.permeability_data,
+ ip_cv.phase_transition_data,
+ ip_dd.dS_L_dp_cap,
+ ip_cv.viscosity_data,
+ ip_dd.dfC_4_LCpC);
+
+ models.fC_4_MCpG_model.dEval(ip_cv.biot_data,
+ current_state.constituent_density_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ ip_dd.porosity_d_data,
+ current_state.S_L_data,
+ ip_cv.beta_p_SR,
+ ip_dd.dfC_4_MCpG);
+
+ models.fC_4_MCT_model.dEval(ip_cv.biot_data,
+ current_state.constituent_density_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ ip_dd.porosity_d_data,
+ current_state.S_L_data,
+ ip_cv.s_therm_exp_data,
+ ip_dd.dfC_4_MCT);
+
+ models.fC_4_MCu_model.dEval(ip_cv.biot_data,
+ ip_cv.phase_transition_data,
+ current_state.S_L_data,
+ ip_dd.dfC_4_MCu);
+
+ if (!this->process_data_.apply_mass_lumping)
{
- double const rho_W_GR_dot =
- (current_state.constituent_density_data.rho_W_GR -
- prev_state.constituent_density_data->rho_W_GR) /
- dt;
- double const rho_W_LR_dot =
- (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt;
-
- ip_cv.dfW_3a_dp_GR =
- /*(ds_G_dp_GR = 0) * rho_W_GR_dot +*/ s_G * c.drho_W_GR_dp_GR /
- dt +
- /*(ds_L_dp_GR = 0) * rho_W_LR_dot +*/ s_L * c.drho_W_LR_dp_GR /
- dt;
- ip_cv.dfW_3a_dp_cap = ds_G_dp_cap * rho_W_GR_dot +
- s_G * c.drho_W_GR_dp_cap / dt +
- ip_cv.dS_L_dp_cap() * rho_W_LR_dot -
- s_L * c.drho_W_LR_dp_LR / dt;
- ip_cv.dfW_3a_dT =
- s_G * c.drho_W_GR_dT / dt + s_L * c.drho_W_LR_dT / dt;
+ models.fW_2_model.dEval(ip_cv.biot_data,
+ pCap_data,
+ current_state.constituent_density_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ ip_dd.porosity_d_data,
+ current_state.rho_W_LR,
+ current_state.S_L_data,
+ ip_dd.dS_L_dp_cap,
+ ip_cv.beta_p_SR,
+ ip_dd.dfW_2);
}
- ip_cv.dfW_4_LWpG_a_dp_GR = c.drho_W_GR_dp_GR * k_over_mu_G
- // + rhoWGR * (dk_over_mu_G_dp_GR = 0)
- + c.drho_W_LR_dp_GR * k_over_mu_L
- // + rhoWLR * (dk_over_mu_L_dp_GR = 0)
- ;
- ip_cv.dfW_4_LWpG_a_dp_cap =
- c.drho_W_GR_dp_cap * k_over_mu_G +
- current_state.constituent_density_data.rho_W_GR *
- ip_cv.dk_over_mu_G_dp_cap +
- -c.drho_W_LR_dp_LR * k_over_mu_L +
- current_state.rho_W_LR() * ip_cv.dk_over_mu_L_dp_cap;
-
- ip_cv.dfW_4_LWpG_a_dT =
- c.drho_W_GR_dT * k_over_mu_G
- //+ rhoWGR * (dk_over_mu_G_dT != 0 TODO for mu_G(T))
- + c.drho_W_LR_dT * k_over_mu_L
- //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_G(T))
- ;
-
- // TODO (naumov) for dxmW*/d* != 0
- ip_cv.dfW_4_LWpG_d_dp_GR =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
- ip_cv.dfW_4_LWpG_d_dp_cap =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
- ip_cv.dfW_4_LWpG_d_dT =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
-
- ip_cv.dfW_4_LWpC_a_dp_GR = c.drho_W_LR_dp_GR * k_over_mu_L
- //+ rhoWLR * (dk_over_mu_L_dp_GR = 0)
- ;
- ip_cv.dfW_4_LWpC_a_dp_cap =
- -c.drho_W_LR_dp_LR * k_over_mu_L +
- current_state.rho_W_LR() * ip_cv.dk_over_mu_L_dp_cap;
- ip_cv.dfW_4_LWpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
- //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T))
- ;
-
- // TODO (naumov) for dxmW*/d* != 0
- ip_cv.dfW_4_LWpC_d_dp_GR =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
- ip_cv.dfW_4_LWpC_d_dp_cap =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
- ip_cv.dfW_4_LWpC_d_dT =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
-
- ip_cv.dfC_4_LCpC_a_dp_GR = c.drho_C_LR_dp_GR * k_over_mu_L
- //+ rhoCLR * (dk_over_mu_L_dp_GR = 0)
- ;
- ip_cv.dfC_4_LCpC_a_dp_cap =
- -c.drho_C_LR_dp_LR * k_over_mu_L +
- current_state.constituent_density_data.rho_C_LR *
- ip_cv.dk_over_mu_L_dp_cap;
- ip_cv.dfC_4_LCpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
- //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T))
- ;
-
- // TODO (naumov) for dxmW*/d* != 0
- ip_cv.dfC_4_LCpC_d_dp_GR =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
- ip_cv.dfC_4_LCpC_d_dp_cap =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
- ip_cv.dfC_4_LCpC_d_dT =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
+ models.fW_3a_model.dEval(dt,
+ current_state.constituent_density_data,
+ ip_cv.phase_transition_data,
+ prev_state.constituent_density_data,
+ prev_state.rho_W_LR,
+ current_state.rho_W_LR,
+ current_state.S_L_data,
+ ip_dd.dS_L_dp_cap,
+ ip_dd.dfW_3a);
+
+ models.fW_4_LWpG_model.dEval(current_state.constituent_density_data,
+ ip_out.permeability_data,
+ ip_cv.phase_transition_data,
+ current_state.rho_W_LR,
+ ip_dd.dS_L_dp_cap,
+ ip_cv.viscosity_data,
+ ip_dd.dfW_4_LWpG);
+
+ models.fW_4_LWpC_model.dEval(ip_cv.advection_data,
+ ip_out.fluid_density_data,
+ ip_out.permeability_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ current_state.rho_W_LR,
+ current_state.S_L_data,
+ ip_dd.dS_L_dp_cap,
+ ip_cv.viscosity_data,
+ ip_dd.dfW_4_LWpC);
+
+ models.fT_1_model.dEval(
+ dt, ip_dd.effective_volumetric_internal_energy_d_data, ip_dd.dfT_1);
+
+ models.fT_2_model.dEval(
+ ip_out.darcy_velocity_data,
+ ip_out.fluid_density_data,
+ ip_out.fluid_enthalpy_data,
+ ip_out.permeability_data,
+ ip_cv.phase_transition_data,
+ ConstitutiveRelations::SpecificBodyForceData<DisplacementDim>{
+ this->process_data_.specific_body_force},
+ ip_cv.viscosity_data,
+ ip_dd.dfT_2);
+
+ models.fu_1_KuT_model.dEval(ip_cd.s_mech_data, ip_cv.s_therm_exp_data,
+ ip_dd.dfu_1_KuT);
+
+ models.fu_2_KupC_model.dEval(ip_cv.biot_data,
+ ip_cv.chi_S_L,
+ pCap_data,
+ ip_dd.dS_L_dp_cap,
+ ip_dd.dfu_2_KupC);
}
- return {ip_constitutive_data, ip_constitutive_variables};
+ return ip_d_data;
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
@@ -1058,10 +1068,6 @@ void TH2MLocalAssembler<
auto const& w = ip.integration_weight;
- auto const& m = Invariants::identity2;
-
- auto const mT = m.transpose().eval();
-
auto const x_coord =
NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement>(
@@ -1073,299 +1079,117 @@ void TH2MLocalAssembler<
typename BMatricesType::BMatrixType>(
gradNu, Nu, x_coord, this->is_axially_symmetric_);
- auto const BuT = Bu.transpose().eval();
+ auto const NTN = (Np.transpose() * Np).eval();
+ auto const BTI2N = (Bu.transpose() * Invariants::identity2 * Np).eval();
double const pCap = Np.dot(capillary_pressure);
double const pCap_prev = Np.dot(capillary_pressure_prev);
- double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
-
- auto const I =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
-
- double const sD_G =
- ip_cv.phase_transition_data.diffusion_coefficient_vapour;
- double const sD_L =
- ip_cv.phase_transition_data.diffusion_coefficient_solute;
-
- GlobalDimMatrixType const D_C_G = sD_G * I;
- GlobalDimMatrixType const D_W_G = sD_G * I;
- GlobalDimMatrixType const D_C_L = sD_L * I;
- GlobalDimMatrixType const D_W_L = sD_L * I;
-
auto const s_L = current_state.S_L_data.S_L;
- auto const s_G = 1. - s_L;
auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
- auto& alpha_B = ip_cv.biot_data();
- auto& beta_p_SR = ip_cv.beta_p_SR();
-
auto const& b = this->process_data_.specific_body_force;
- // volume fraction
- auto const phi_G = s_G * ip_out.porosity_data.phi;
- auto const phi_L = s_L * ip_out.porosity_data.phi;
- auto const phi_S = 1. - ip_out.porosity_data.phi;
-
- auto const rhoGR = ip_out.fluid_density_data.rho_GR;
- auto const rhoLR = ip_out.fluid_density_data.rho_LR;
-
- // effective density
- auto const rho = phi_G * rhoGR + phi_L * rhoLR +
- phi_S * ip_out.solid_density_data.rho_SR;
-
- // abbreviations
- const double rho_C_FR =
- s_G * current_state.constituent_density_data.rho_C_GR +
- s_L * current_state.constituent_density_data.rho_C_LR;
- const double rho_W_FR =
- s_G * current_state.constituent_density_data.rho_W_GR +
- s_L * current_state.rho_W_LR();
-
- auto const rho_C_GR_dot =
- (current_state.constituent_density_data.rho_C_GR -
- prev_state.constituent_density_data->rho_C_GR) /
- dt;
- auto const rho_C_LR_dot =
- (current_state.constituent_density_data.rho_C_LR -
- prev_state.constituent_density_data->rho_C_LR) /
- dt;
- auto const rho_W_GR_dot =
- (current_state.constituent_density_data.rho_W_GR -
- prev_state.constituent_density_data->rho_W_GR) /
- dt;
- auto const rho_W_LR_dot =
- (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt;
-
- GlobalDimMatrixType const k_over_mu_G =
- ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
- ip_cv.viscosity_data.mu_GR;
- GlobalDimMatrixType const k_over_mu_L =
- ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
- ip_cv.viscosity_data.mu_LR;
-
// ---------------------------------------------------------------------
// C-component equation
// ---------------------------------------------------------------------
- MCpG.noalias() += NpT * rho_C_FR *
- (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
- Np * w;
- MCpC.noalias() -= NpT * rho_C_FR *
- (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
- s_L * Np * w;
+ MCpG.noalias() += NTN * (ip_cv.fC_4_MCpG.m * w);
+ MCpC.noalias() += NTN * (ip_cv.fC_4_MCpC.m * w);
if (this->process_data_.apply_mass_lumping)
{
if (pCap - pCap_prev != 0.) // avoid division by Zero
{
MCpC.noalias() +=
- NpT *
- (ip_out.porosity_data.phi *
- (current_state.constituent_density_data.rho_C_LR -
- current_state.constituent_density_data.rho_C_GR) -
- rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
- beta_p_SR) *
- s_L_dot * dt / (pCap - pCap_prev) * Np * w;
+ NTN * (ip_cv.fC_4_MCpC.ml / (pCap - pCap_prev) * w);
}
}
- MCT.noalias() -= NpT * rho_C_FR * (alpha_B - ip_out.porosity_data.phi) *
- beta_T_SR * Np * w;
- MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
-
- using DisplacementDimMatrix =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
-
- DisplacementDimMatrix const advection_C_G =
- current_state.constituent_density_data.rho_C_GR * k_over_mu_G;
- DisplacementDimMatrix const advection_C_L =
- current_state.constituent_density_data.rho_C_LR * k_over_mu_L;
-
- DisplacementDimMatrix const diffusion_CGpGR =
- -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpGR;
- DisplacementDimMatrix const diffusion_CLpGR =
- -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpGR;
-
- DisplacementDimMatrix const diffusion_CGpCap =
- -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpCap;
- DisplacementDimMatrix const diffusion_CLpCap =
- -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpCap;
-
- DisplacementDimMatrix const diffusion_CGT =
- -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dT;
- DisplacementDimMatrix const diffusion_CLT =
- -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dT;
+ MCT.noalias() += NTN * (ip_cv.fC_4_MCT.m * w);
+ MCu.noalias() += BTI2N.transpose() * (ip_cv.fC_4_MCu.m * w);
- DisplacementDimMatrix const advection_C = advection_C_G + advection_C_L;
- DisplacementDimMatrix const diffusion_C_pGR =
- diffusion_CGpGR + diffusion_CLpGR;
- DisplacementDimMatrix const diffusion_C_pCap =
- diffusion_CGpCap + diffusion_CLpCap;
+ LCpG.noalias() += gradNpT * ip_cv.fC_4_LCpG.L * gradNp * w;
- DisplacementDimMatrix const diffusion_C_T =
- diffusion_CGT + diffusion_CLT;
+ LCpC.noalias() += gradNpT * ip_cv.fC_4_LCpC.L * gradNp * w;
- LCpG.noalias() +=
- gradNpT * (advection_C + diffusion_C_pGR) * gradNp * w;
+ LCT.noalias() += gradNpT * ip_cv.fC_4_LCT.L * gradNp * w;
- LCpC.noalias() +=
- gradNpT * (diffusion_C_pCap - advection_C_L) * gradNp * w;
-
- LCT.noalias() += gradNpT * (diffusion_C_T)*gradNp * w;
-
- fC.noalias() +=
- gradNpT * (advection_C_G * rhoGR + advection_C_L * rhoLR) * b * w;
+ fC.noalias() += gradNpT * ip_cv.fC_1.A * b * w;
if (!this->process_data_.apply_mass_lumping)
{
- fC.noalias() -=
- NpT *
- (ip_out.porosity_data.phi *
- (current_state.constituent_density_data.rho_C_LR -
- current_state.constituent_density_data.rho_C_GR) -
- rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
- beta_p_SR) *
- s_L_dot * w;
+ fC.noalias() -= NpT * (ip_cv.fC_2a.a * s_L_dot * w);
}
// fC_III
- fC.noalias() -= NpT * ip_out.porosity_data.phi *
- (s_G * rho_C_GR_dot + s_L * rho_C_LR_dot) * w;
+ fC.noalias() -= NpT * (ip_out.porosity_data.phi * ip_cv.fC_3a.a * w);
// ---------------------------------------------------------------------
// W-component equation
// ---------------------------------------------------------------------
- MWpG.noalias() += NpT * rho_W_FR *
- (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
- Np * w;
- MWpC.noalias() -= NpT * rho_W_FR *
- (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
- s_L * Np * w;
+ MWpG.noalias() += NTN * (ip_cv.fW_4_MWpG.m * w);
+ MWpC.noalias() += NTN * (ip_cv.fW_4_MWpC.m * w);
if (this->process_data_.apply_mass_lumping)
{
if (pCap - pCap_prev != 0.) // avoid division by Zero
{
MWpC.noalias() +=
- NpT *
- (ip_out.porosity_data.phi *
- (current_state.rho_W_LR() -
- current_state.constituent_density_data.rho_W_GR) -
- rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
- beta_p_SR) *
- s_L_dot * dt / (pCap - pCap_prev) * Np * w;
+ NTN * (ip_cv.fW_4_MWpC.ml / (pCap - pCap_prev) * w);
}
}
- MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
- beta_T_SR * Np * w;
-
- MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
-
- DisplacementDimMatrix const advection_W_G =
- current_state.constituent_density_data.rho_W_GR * k_over_mu_G;
- DisplacementDimMatrix const advection_W_L =
- current_state.rho_W_LR() * k_over_mu_L;
-
- DisplacementDimMatrix const diffusion_WGpGR =
- phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpGR;
- DisplacementDimMatrix const diffusion_WLpGR =
- phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpGR;
-
- DisplacementDimMatrix const diffusion_WGpCap =
- phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpCap;
- DisplacementDimMatrix const diffusion_WLpCap =
- phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpCap;
-
- DisplacementDimMatrix const diffusion_WGT =
- phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dT;
- DisplacementDimMatrix const diffusion_WLT =
- phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dT;
-
- DisplacementDimMatrix const advection_W = advection_W_G + advection_W_L;
- DisplacementDimMatrix const diffusion_W_pGR =
- diffusion_WGpGR + diffusion_WLpGR;
- DisplacementDimMatrix const diffusion_W_pCap =
- diffusion_WGpCap + diffusion_WLpCap;
+ MWT.noalias() += NTN * (ip_cv.fW_4_MWT.m * w);
- DisplacementDimMatrix const diffusion_W_T =
- diffusion_WGT + diffusion_WLT;
+ MWu.noalias() += BTI2N.transpose() * (ip_cv.fW_4_MWu.m * w);
- LWpG.noalias() +=
- gradNpT * (advection_W + diffusion_W_pGR) * gradNp * w;
+ LWpG.noalias() += gradNpT * ip_cv.fW_4_LWpG.L * gradNp * w;
- LWpC.noalias() +=
- gradNpT * (diffusion_W_pCap - advection_W_L) * gradNp * w;
+ LWpC.noalias() += gradNpT * ip_cv.fW_4_LWpC.L * gradNp * w;
- LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
+ LWT.noalias() += gradNpT * ip_cv.fW_4_LWT.L * gradNp * w;
- fW.noalias() +=
- gradNpT * (advection_W_G * rhoGR + advection_W_L * rhoLR) * b * w;
+ fW.noalias() += gradNpT * ip_cv.fW_1.A * b * w;
if (!this->process_data_.apply_mass_lumping)
{
- fW.noalias() -=
- NpT *
- (ip_out.porosity_data.phi *
- (current_state.rho_W_LR() -
- current_state.constituent_density_data.rho_W_GR) -
- rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
- beta_p_SR) *
- s_L_dot * w;
+ fW.noalias() -= NpT * (ip_cv.fW_2.a * s_L_dot * w);
}
- fW.noalias() -= NpT * ip_out.porosity_data.phi *
- (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
+ fW.noalias() -= NpT * (ip_out.porosity_data.phi * ip_cv.fW_3a.a * w);
// ---------------------------------------------------------------------
// - temperature equation
// ---------------------------------------------------------------------
- MTu.noalias() += NTT *
- ip_cv.effective_volumetric_enthalpy_data.rho_h_eff *
- mT * Bu * w;
+ MTu.noalias() +=
+ BTI2N.transpose() *
+ (ip_cv.effective_volumetric_enthalpy_data.rho_h_eff * w);
KTT.noalias() +=
gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
- auto const rho_u_eff_dot = (current_state.internal_energy_data() -
- **prev_state.internal_energy_data) /
- dt;
- fT.noalias() -= NTT * rho_u_eff_dot * w;
-
- fT.noalias() += gradNTT *
- (rhoGR * ip_out.enthalpy_data.h_G *
- ip_out.darcy_velocity_data.w_GS +
- rhoLR * ip_out.enthalpy_data.h_L *
- ip_out.darcy_velocity_data.w_LS) *
- w;
-
- fT.noalias() += gradNTT *
- (current_state.constituent_density_data.rho_C_GR *
- ip_cv.phase_transition_data.hCG *
- ip_out.diffusion_velocity_data.d_CG +
- current_state.constituent_density_data.rho_W_GR *
- ip_cv.phase_transition_data.hWG *
- ip_out.diffusion_velocity_data.d_WG) *
- w;
-
- fT.noalias() += NTT *
- (rhoGR * ip_out.darcy_velocity_data.w_GS.transpose() +
- rhoLR * ip_out.darcy_velocity_data.w_LS.transpose()) *
- b * w;
+ fT.noalias() -= NTT * (ip_cv.fT_1.m * w);
+
+ fT.noalias() += gradNTT * ip_cv.fT_2.A * w;
+
+ fT.noalias() += gradNTT * ip_cv.fT_3.gradN * w;
+
+ fT.noalias() += NTT * (ip_cv.fT_3.N * w);
// ---------------------------------------------------------------------
// - displacement equation
// ---------------------------------------------------------------------
- KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
+ KUpG.noalias() -= BTI2N * (ip_cv.biot_data() * w);
- KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w;
+ KUpC.noalias() += BTI2N * (ip_cv.fu_2_KupC.m * w);
- fU.noalias() -= (BuT * current_state.eff_stress_data.sigma -
- N_u_op(Nu).transpose() * rho * b) *
- w;
+ fU.noalias() -=
+ (Bu.transpose() * current_state.eff_stress_data.sigma -
+ N_u_op(Nu).transpose() * ip_cv.volumetric_body_force()) *
+ w;
if (this->process_data_.apply_mass_lumping)
{
@@ -1515,10 +1339,17 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
local_x_prev.size()),
t, dt, models);
+ auto const ip_d_data = updateConstitutiveVariablesDerivatives(
+ Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
+ Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
+ local_x_prev.size()),
+ t, dt, ip_constitutive_data, ip_constitutive_variables, models);
+
for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
{
auto& ip = _ip_data[int_point];
auto& ip_cd = ip_constitutive_data[int_point];
+ auto& ip_dd = ip_d_data[int_point];
auto& ip_cv = ip_constitutive_variables[int_point];
auto& ip_out = this->output_data_[int_point];
auto& current_state = this->current_states_[int_point];
@@ -1546,9 +1377,6 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
auto const& w = ip.integration_weight;
- auto const& m = Invariants::identity2;
- auto const mT = m.transpose().eval();
-
auto const x_coord =
NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement>(
@@ -1560,7 +1388,8 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
typename BMatricesType::BMatrixType>(
gradNu, Nu, x_coord, this->is_axially_symmetric_);
- auto const BuT = Bu.transpose().eval();
+ auto const NTN = (Np.transpose() * Np).eval();
+ auto const BTI2N = (Bu.transpose() * Invariants::identity2 * Np).eval();
double const div_u_dot =
Invariants::trace(Bu * (displacement - displacement_prev) / dt);
@@ -1576,402 +1405,232 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
double const pGR_prev = Np.dot(gas_pressure_prev);
double const pCap_prev = Np.dot(capillary_pressure_prev);
double const T_prev = NT.dot(temperature_prev);
- double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
- auto const I =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
-
- double const sD_G =
- ip_cv.phase_transition_data.diffusion_coefficient_vapour;
- double const sD_L =
- ip_cv.phase_transition_data.diffusion_coefficient_solute;
-
- GlobalDimMatrixType const D_C_G = sD_G * I;
- GlobalDimMatrixType const D_W_G = sD_G * I;
- GlobalDimMatrixType const D_C_L = sD_L * I;
- GlobalDimMatrixType const D_W_L = sD_L * I;
-
- auto& s_L = current_state.S_L_data.S_L;
- auto const s_G = 1. - s_L;
+ auto const& s_L = current_state.S_L_data.S_L;
auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
- auto const alpha_B = ip_cv.biot_data();
- auto const beta_p_SR = ip_cv.beta_p_SR();
-
auto const& b = this->process_data_.specific_body_force;
- // volume fraction
- auto const phi_G = s_G * ip_out.porosity_data.phi;
- auto const phi_L = s_L * ip_out.porosity_data.phi;
- auto const phi_S = 1. - ip_out.porosity_data.phi;
-
- auto const rhoGR = ip_out.fluid_density_data.rho_GR;
- auto const rhoLR = ip_out.fluid_density_data.rho_LR;
-
- // effective density
- auto const rho = phi_G * rhoGR + phi_L * rhoLR +
- phi_S * ip_out.solid_density_data.rho_SR;
-
- // abbreviations
- const double rho_C_FR =
- s_G * current_state.constituent_density_data.rho_C_GR +
- s_L * current_state.constituent_density_data.rho_C_LR;
- const double rho_W_FR =
- s_G * current_state.constituent_density_data.rho_W_GR +
- s_L * current_state.rho_W_LR();
-
- auto const rho_C_GR_dot =
- (current_state.constituent_density_data.rho_C_GR -
- prev_state.constituent_density_data->rho_C_GR) /
- dt;
- auto const rho_C_LR_dot =
- (current_state.constituent_density_data.rho_C_LR -
- prev_state.constituent_density_data->rho_C_LR) /
- dt;
- auto const rho_W_GR_dot =
- (current_state.constituent_density_data.rho_W_GR -
- prev_state.constituent_density_data->rho_W_GR) /
- dt;
- auto const rho_W_LR_dot =
- (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt;
-
- GlobalDimMatrixType const k_over_mu_G =
- ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
- ip_cv.viscosity_data.mu_GR;
- GlobalDimMatrixType const k_over_mu_L =
- ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
- ip_cv.viscosity_data.mu_LR;
-
// ---------------------------------------------------------------------
// C-component equation
// ---------------------------------------------------------------------
- MCpG.noalias() += NpT * rho_C_FR *
- (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
- Np * w;
- MCpC.noalias() -= NpT * rho_C_FR *
- (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
- s_L * Np * w;
+ MCpG.noalias() += NTN * (ip_cv.fC_4_MCpG.m * w);
+ MCpC.noalias() += NTN * (ip_cv.fC_4_MCpC.m * w);
if (this->process_data_.apply_mass_lumping)
{
if (pCap - pCap_prev != 0.) // avoid division by Zero
{
MCpC.noalias() +=
- NpT *
- (ip_out.porosity_data.phi *
- (current_state.constituent_density_data.rho_C_LR -
- current_state.constituent_density_data.rho_C_GR) -
- rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
- beta_p_SR) *
- s_L_dot * dt / (pCap - pCap_prev) * Np * w;
+ NTN * (ip_cv.fC_4_MCpC.ml / (pCap - pCap_prev) * w);
}
}
- MCT.noalias() -= NpT * rho_C_FR * (alpha_B - ip_out.porosity_data.phi) *
- beta_T_SR * Np * w;
+ MCT.noalias() += NTN * (ip_cv.fC_4_MCT.m * w);
// d (fC_4_MCT * T_dot)/d T
local_Jac
.template block<C_size, temperature_size>(C_index,
temperature_index)
- .noalias() += NpT * ip_cv.dfC_4_MCT_dT * (T - T_prev) / dt * NT * w;
+ .noalias() += NTN * (ip_dd.dfC_4_MCT.dT * (T - T_prev) / dt * w);
- MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
+ MCu.noalias() += BTI2N.transpose() * (ip_cv.fC_4_MCu.m * w);
// d (fC_4_MCu * u_dot)/d T
local_Jac
.template block<C_size, temperature_size>(C_index,
temperature_index)
- .noalias() += NpT * ip_cv.dfC_4_MCu_dT * div_u_dot * NT * w;
-
- GlobalDimMatrixType const advection_C_G =
- current_state.constituent_density_data.rho_C_GR * k_over_mu_G;
- GlobalDimMatrixType const advection_C_L =
- current_state.constituent_density_data.rho_C_LR * k_over_mu_L;
- GlobalDimMatrixType const diffusion_C_G_p =
- -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpGR;
- GlobalDimMatrixType const diffusion_C_L_p =
- -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpLR;
- GlobalDimMatrixType const diffusion_C_G_T =
- -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dT;
- GlobalDimMatrixType const diffusion_C_L_T =
- -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dT;
-
- GlobalDimMatrixType const advection_C = advection_C_G + advection_C_L;
- GlobalDimMatrixType const diffusion_C_p =
- diffusion_C_G_p + diffusion_C_L_p;
- GlobalDimMatrixType const diffusion_C_T =
- diffusion_C_G_T + diffusion_C_L_T;
-
- LCpG.noalias() += gradNpT * (advection_C + diffusion_C_p) * gradNp * w;
+ .noalias() += NTN * (ip_dd.dfC_4_MCu.dT * div_u_dot * w);
+
+ LCpG.noalias() += gradNpT * ip_cv.fC_4_LCpG.L * gradNp * w;
// d (fC_4_LCpG * grad p_GR)/d p_GR
local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
- gradNpT *
- (ip_cv.dadvection_C_dp_GR
- // + ip_cv.ddiffusion_C_p_dp_GR TODO (naumov)
- ) *
- gradpGR * Np * w;
+ gradNpT * ip_dd.dfC_4_LCpG.dp_GR * gradpGR * Np * w;
// d (fC_4_LCpG * grad p_GR)/d p_cap
local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
- gradNpT *
- (ip_cv.dadvection_C_dp_cap
- // + ip_cv.ddiffusion_C_p_dp_GR TODO (naumov)
- ) *
- gradpGR * Np * w;
+ gradNpT * ip_dd.dfC_4_LCpG.dp_cap * gradpGR * Np * w;
// d (fC_4_LCpG * grad p_GR)/d T
local_Jac
.template block<C_size, temperature_size>(C_index,
temperature_index)
- .noalias() += gradNpT * ip_cv.dfC_4_LCpG_dT * gradpGR * NT * w;
+ .noalias() += gradNpT * ip_dd.dfC_4_LCpG.dT * gradpGR * NT * w;
// d (fC_4_MCpG * p_GR_dot)/d p_GR
local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
- NpT * ip_cv.dfC_4_MCpG_dp_GR * (pGR - pGR_prev) / dt * Np * w;
+ NTN * (ip_dd.dfC_4_MCpG.dp_GR * (pGR - pGR_prev) / dt * w);
// d (fC_4_MCpG * p_GR_dot)/d T
local_Jac
.template block<C_size, temperature_size>(C_index,
temperature_index)
.noalias() +=
- NpT * ip_cv.dfC_4_MCpG_dT * (pGR - pGR_prev) / dt * NT * w;
+ NTN * (ip_dd.dfC_4_MCpG.dT * (pGR - pGR_prev) / dt * w);
- LCpC.noalias() -=
- gradNpT * (advection_C_L + diffusion_C_L_p) * gradNp * w;
+ LCpC.noalias() -= gradNpT * ip_cv.fC_4_LCpC.L * gradNp * w;
/* TODO (naumov) This part is not tested by any of the current ctests.
// d (fC_4_LCpC * grad p_cap)/d p_GR
local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
- gradNpT *
- (ip_cv.dfC_4_LCpC_a_dp_GR
- // + ip_cv.dfC_4_LCpC_d_dp_GR TODO (naumov)
- ) *
- gradpCap * Np * w;
+ gradNpT * ip_dd.dfC_4_LCpC.dp_GR * gradpCap * Np * w;
// d (fC_4_LCpC * grad p_cap)/d p_cap
local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
- gradNpT *
- (ip_cv.dfC_4_LCpC_a_dp_cap
- // + ip_cv.dfC_4_LCpC_d_dp_cap TODO (naumov)
- ) *
- gradpCap * Np * w;
+ gradNpT * ip_dd.dfC_4_LCpC.dp_cap * gradpCap * Np * w;
local_Jac
.template block<C_size, temperature_size>(C_index,
temperature_index)
- .noalias() += gradNpT *
- (ip_cv.dfC_4_LCpC_a_dT
- // + ip_cv.dfC_4_LCpC_d_dT TODO (naumov)
- ) *
- gradpCap * Np * w;
+ .noalias() += gradNpT * ip_dd.dfC_4_LCpC.dT * gradpCap * Np * w;
*/
- LCT.noalias() += gradNpT * diffusion_C_T * gradNp * w;
+ LCT.noalias() += gradNpT * ip_cv.fC_4_LCT.L * gradNp * w;
// fC_1
- fC.noalias() +=
- gradNpT * (advection_C_G * rhoGR + advection_C_L * rhoLR) * b * w;
+ fC.noalias() += gradNpT * ip_cv.fC_1.A * b * w;
if (!this->process_data_.apply_mass_lumping)
{
// fC_2 = \int a * s_L_dot
- auto const a =
- ip_out.porosity_data.phi *
- (current_state.constituent_density_data.rho_C_LR -
- current_state.constituent_density_data.rho_C_GR) -
- rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
- beta_p_SR;
- fC.noalias() -= NpT * a * s_L_dot * w;
+ fC.noalias() -= NpT * (ip_cv.fC_2a.a * s_L_dot * w);
local_Jac.template block<C_size, C_size>(C_index, C_index)
.noalias() +=
- NpT *
- (ip_cv.dfC_2a_dp_GR * s_L_dot /*- a * (ds_L_dp_GR = 0) / dt*/) *
- Np * w;
+ NTN * ((ip_dd.dfC_2a.dp_GR * s_L_dot
+ /*- ip_cv.fC_2a.a * (ds_L_dp_GR = 0) / dt*/) *
+ w);
local_Jac.template block<C_size, W_size>(C_index, W_index)
.noalias() +=
- NpT *
- (ip_cv.dfC_2a_dp_cap * s_L_dot + a * ip_cv.dS_L_dp_cap() / dt) *
- Np * w;
+ NTN * ((ip_dd.dfC_2a.dp_cap * s_L_dot +
+ ip_cv.fC_2a.a * ip_dd.dS_L_dp_cap() / dt) *
+ w);
local_Jac
.template block<C_size, temperature_size>(C_index,
temperature_index)
- .noalias() += NpT * ip_cv.dfC_2a_dT * s_L_dot * NT * w;
+ .noalias() += NTN * (ip_dd.dfC_2a.dT * s_L_dot * w);
}
{
// fC_3 = \int phi * a
- double const a = s_G * rho_C_GR_dot + s_L * rho_C_LR_dot;
- fC.noalias() -= NpT * ip_out.porosity_data.phi * a * w;
+ fC.noalias() -=
+ NpT * (ip_out.porosity_data.phi * ip_cv.fC_3a.a * w);
local_Jac.template block<C_size, C_size>(C_index, C_index)
.noalias() +=
- NpT * ip_out.porosity_data.phi * ip_cv.dfC_3a_dp_GR * Np * w;
+ NTN * (ip_out.porosity_data.phi * ip_dd.dfC_3a.dp_GR * w);
local_Jac.template block<C_size, W_size>(C_index, W_index)
.noalias() +=
- NpT * ip_out.porosity_data.phi * ip_cv.dfC_3a_dp_cap * Np * w;
+ NTN * (ip_out.porosity_data.phi * ip_dd.dfC_3a.dp_cap * w);
local_Jac
.template block<C_size, temperature_size>(C_index,
temperature_index)
- .noalias() += NpT *
- (ip_cv.porosity_d_data.dphi_dT * a +
- ip_out.porosity_data.phi * ip_cv.dfC_3a_dT) *
- NT * w;
+ .noalias() +=
+ NTN * ((ip_dd.porosity_d_data.dphi_dT * ip_cv.fC_3a.a +
+ ip_out.porosity_data.phi * ip_dd.dfC_3a.dT) *
+ w);
}
// ---------------------------------------------------------------------
// W-component equation
// ---------------------------------------------------------------------
- MWpG.noalias() += NpT * rho_W_FR *
- (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
- Np * w;
- MWpC.noalias() -= NpT * rho_W_FR *
- (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
- s_L * Np * w;
+ MWpG.noalias() += NTN * (ip_cv.fW_4_MWpG.m * w);
+ MWpC.noalias() += NTN * (ip_cv.fW_4_MWpC.m * w);
if (this->process_data_.apply_mass_lumping)
{
if (pCap - pCap_prev != 0.) // avoid division by Zero
{
MWpC.noalias() +=
- NpT *
- (ip_out.porosity_data.phi *
- (current_state.rho_W_LR() -
- current_state.constituent_density_data.rho_W_GR) -
- rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
- beta_p_SR) *
- s_L_dot * dt / (pCap - pCap_prev) * Np * w;
+ NTN * (ip_cv.fW_4_MWpC.ml / (pCap - pCap_prev) * w);
}
}
- MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
- beta_T_SR * Np * w;
-
- MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
+ MWT.noalias() += NTN * (ip_cv.fW_4_MWT.m * w);
- GlobalDimMatrixType const advection_W_G =
- current_state.constituent_density_data.rho_W_GR * k_over_mu_G;
- GlobalDimMatrixType const advection_W_L =
- current_state.rho_W_LR() * k_over_mu_L;
- GlobalDimMatrixType const diffusion_W_G_p =
- phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpGR;
- GlobalDimMatrixType const diffusion_W_L_p =
- phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpLR;
- GlobalDimMatrixType const diffusion_W_G_T =
- phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dT;
- GlobalDimMatrixType const diffusion_W_L_T =
- phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dT;
+ MWu.noalias() += BTI2N.transpose() * (ip_cv.fW_4_MWu.m * w);
- GlobalDimMatrixType const advection_W = advection_W_G + advection_W_L;
- GlobalDimMatrixType const diffusion_W_p =
- diffusion_W_G_p + diffusion_W_L_p;
- GlobalDimMatrixType const diffusion_W_T =
- diffusion_W_G_T + diffusion_W_L_T;
-
- LWpG.noalias() += gradNpT * (advection_W + diffusion_W_p) * gradNp * w;
+ LWpG.noalias() += gradNpT * ip_cv.fW_4_LWpG.L * gradNp * w;
// fW_4 LWpG' parts; LWpG = \int grad (a + d) grad
local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
- gradNpT * (ip_cv.dfW_4_LWpG_a_dp_GR + ip_cv.dfW_4_LWpG_d_dp_GR) *
- gradpGR * Np * w;
+ gradNpT * ip_dd.dfW_4_LWpG.dp_GR * gradpGR * Np * w;
local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
- gradNpT * (ip_cv.dfW_4_LWpG_a_dp_cap + ip_cv.dfW_4_LWpG_d_dp_cap) *
- gradpGR * Np * w;
+ gradNpT * ip_dd.dfW_4_LWpG.dp_cap * gradpGR * Np * w;
local_Jac
.template block<W_size, temperature_size>(W_index,
temperature_index)
- .noalias() += gradNpT *
- (ip_cv.dfW_4_LWpG_a_dT + ip_cv.dfW_4_LWpG_d_dT) *
- gradpGR * NT * w;
+ .noalias() += gradNpT * ip_dd.dfW_4_LWpG.dT * gradpGR * NT * w;
- LWpC.noalias() -=
- gradNpT * (advection_W_L + diffusion_W_L_p) * gradNp * w;
+ LWpC.noalias() += gradNpT * ip_cv.fW_4_LWpC.L * gradNp * w;
// fW_4 LWp_cap' parts; LWpC = \int grad (a + d) grad
local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
- gradNpT * (ip_cv.dfW_4_LWpC_a_dp_GR + ip_cv.dfW_4_LWpC_d_dp_GR) *
- gradpCap * Np * w;
+ gradNpT * ip_dd.dfW_4_LWpC.dp_GR * gradpCap * Np * w;
local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
- gradNpT * (ip_cv.dfW_4_LWpC_a_dp_cap + ip_cv.dfW_4_LWpC_d_dp_cap) *
- gradpCap * Np * w;
+ gradNpT * ip_dd.dfW_4_LWpC.dp_cap * gradpCap * Np * w;
local_Jac
.template block<W_size, temperature_size>(W_index,
temperature_index)
- .noalias() -= gradNpT *
- (ip_cv.dfW_4_LWpC_a_dT + ip_cv.dfW_4_LWpC_d_dT) *
- gradpCap * NT * w;
+ .noalias() -= gradNpT * ip_dd.dfW_4_LWpC.dT * gradpCap * NT * w;
- LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
+ LWT.noalias() += gradNpT * ip_cv.fW_4_LWT.L * gradNp * w;
// fW_1
- fW.noalias() +=
- gradNpT * (advection_W_G * rhoGR + advection_W_L * rhoLR) * b * w;
+ fW.noalias() += gradNpT * ip_cv.fW_1.A * b * w;
- // fW_2 = \int (f - g) * s_L_dot
+ // fW_2 = \int a * s_L_dot
if (!this->process_data_.apply_mass_lumping)
{
- double const f = ip_out.porosity_data.phi *
- (current_state.rho_W_LR() -
- current_state.constituent_density_data.rho_W_GR);
- double const g = rho_W_FR * pCap *
- (alpha_B - ip_out.porosity_data.phi) * beta_p_SR;
-
- fW.noalias() -= NpT * (f - g) * s_L_dot * w;
+ fW.noalias() -= NpT * (ip_cv.fW_2.a * s_L_dot * w);
local_Jac.template block<W_size, C_size>(W_index, C_index)
- .noalias() += NpT * (ip_cv.dfW_2a_dp_GR - ip_cv.dfW_2b_dp_GR) *
- s_L_dot * Np * w;
+ .noalias() += NTN * (ip_dd.dfW_2.dp_GR * s_L_dot * w);
// sign negated because of dp_cap = -dp_LR
// TODO (naumov) Had to change the sign to get equal Jacobian WW
// blocks in A2 Test. Where is the error?
local_Jac.template block<W_size, W_size>(W_index, W_index)
- .noalias() +=
- NpT *
- ((ip_cv.dfW_2a_dp_cap - ip_cv.dfW_2b_dp_cap) * s_L_dot +
- (f - g) * ip_cv.dS_L_dp_cap() / dt) *
- Np * w;
+ .noalias() += NTN * ((ip_dd.dfW_2.dp_cap * s_L_dot +
+ ip_cv.fW_2.a * ip_dd.dS_L_dp_cap() / dt) *
+ w);
local_Jac
.template block<W_size, temperature_size>(W_index,
temperature_index)
- .noalias() +=
- NpT * (ip_cv.dfW_2a_dT - ip_cv.dfW_2b_dT) * s_L_dot * Np * w;
+ .noalias() += NTN * (ip_dd.dfW_2.dT * s_L_dot * w);
}
// fW_3 = \int phi * a
- fW.noalias() -= NpT * ip_out.porosity_data.phi *
- (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
+ fW.noalias() -= NpT * (ip_out.porosity_data.phi * ip_cv.fW_3a.a * w);
local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
- NpT * ip_out.porosity_data.phi * ip_cv.dfW_3a_dp_GR * Np * w;
+ NTN * (ip_out.porosity_data.phi * ip_dd.dfW_3a.dp_GR * w);
local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
- NpT * ip_out.porosity_data.phi * ip_cv.dfW_3a_dp_cap * Np * w;
+ NTN * (ip_out.porosity_data.phi * ip_dd.dfW_3a.dp_cap * w);
local_Jac
.template block<W_size, temperature_size>(W_index,
temperature_index)
- .noalias() += NpT *
- (ip_cv.porosity_d_data.dphi_dT *
- (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) +
- ip_out.porosity_data.phi * ip_cv.dfW_3a_dT) *
- NT * w;
+ .noalias() +=
+ NTN * ((ip_dd.porosity_d_data.dphi_dT * ip_cv.fW_3a.a +
+ ip_out.porosity_data.phi * ip_dd.dfW_3a.dT) *
+ w);
// ---------------------------------------------------------------------
// - temperature equation
// ---------------------------------------------------------------------
- MTu.noalias() += NTT *
- ip_cv.effective_volumetric_enthalpy_data.rho_h_eff *
- mT * Bu * w;
+ MTu.noalias() +=
+ BTI2N.transpose() *
+ (ip_cv.effective_volumetric_enthalpy_data.rho_h_eff * w);
// dfT_4/dp_GR
// d (MTu * u_dot)/dp_GR
@@ -1979,8 +1638,8 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
.template block<temperature_size, C_size>(temperature_index,
C_index)
.noalias() +=
- NTT * ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_GR *
- div_u_dot * NT * w;
+ NTN * (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_GR *
+ div_u_dot * w);
// dfT_4/dp_cap
// d (MTu * u_dot)/dp_cap
@@ -1988,8 +1647,9 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
.template block<temperature_size, W_size>(temperature_index,
W_index)
.noalias() -=
- NTT * ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_cap *
- div_u_dot * NT * w;
+ NTN *
+ (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_cap *
+ div_u_dot * w);
// dfT_4/dT
// d (MTu * u_dot)/dT
@@ -1997,8 +1657,8 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
.template block<temperature_size, temperature_size>(
temperature_index, temperature_index)
.noalias() +=
- NTT * ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dT *
- div_u_dot * NT * w;
+ NTN * (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dT *
+ div_u_dot * w);
KTT.noalias() +=
gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
@@ -2022,57 +1682,41 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
.template block<temperature_size, W_size>(temperature_index,
W_index)
.noalias() += gradNTT *
- ip_cv.thermal_conductivity_data.dlambda_dp_cap *
+ ip_dd.thermal_conductivity_d_data.dlambda_dp_cap *
gradT * Np * w;
// d KTT/dT * T
local_Jac
.template block<temperature_size, temperature_size>(
temperature_index, temperature_index)
- .noalias() += gradNTT * ip_cv.thermal_conductivity_data.dlambda_dT *
- gradT * NT * w;
+ .noalias() += gradNTT *
+ ip_dd.thermal_conductivity_d_data.dlambda_dT * gradT *
+ NT * w;
// fT_1
- auto const rho_u_eff_dot = (current_state.internal_energy_data() -
- **prev_state.internal_energy_data) /
- dt;
- fT.noalias() -= NTT * rho_u_eff_dot * w;
+ fT.noalias() -= NTT * (ip_cv.fT_1.m * w);
// dfT_1/dp_GR
local_Jac
.template block<temperature_size, C_size>(temperature_index,
C_index)
- .noalias() += NpT * Np *
- (ip_cv.effective_volumetric_internal_energy_d_data
- .drho_u_eff_dp_GR /
- dt * w);
+ .noalias() += NTN * (ip_dd.dfT_1.dp_GR * w);
// dfT_1/dp_cap
local_Jac
.template block<temperature_size, W_size>(temperature_index,
W_index)
- .noalias() += NpT * Np *
- (ip_cv.effective_volumetric_internal_energy_d_data
- .drho_u_eff_dp_cap /
- dt * w);
+ .noalias() += NTN * (ip_dd.dfT_1.dp_cap * w);
// dfT_1/dT
// MTT
local_Jac
.template block<temperature_size, temperature_size>(
temperature_index, temperature_index)
- .noalias() +=
- NTT * NT *
- (ip_cv.effective_volumetric_internal_energy_d_data.drho_u_eff_dT /
- dt * w);
+ .noalias() += NTN * (ip_dd.dfT_1.dT * w);
// fT_2
- fT.noalias() += gradNTT *
- (rhoGR * ip_out.enthalpy_data.h_G *
- ip_out.darcy_velocity_data.w_GS +
- rhoLR * ip_out.enthalpy_data.h_L *
- ip_out.darcy_velocity_data.w_LS) *
- w;
+ fT.noalias() += gradNTT * ip_cv.fT_2.A * w;
// dfT_2/dp_GR
local_Jac
@@ -2080,9 +1724,9 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
C_index)
.noalias() -=
// dfT_2/dp_GR first part
- gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_Npart * Np * w +
+ gradNTT * ip_dd.dfT_2.dp_GR_Npart * Np * w +
// dfT_2/dp_GR second part
- gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart * gradNp * w;
+ gradNTT * ip_dd.dfT_2.dp_GR_gradNpart * gradNp * w;
// dfT_2/dp_cap
local_Jac
@@ -2090,69 +1734,56 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
W_index)
.noalias() -=
// first part of dfT_2/dp_cap
- gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_Npart) * Np * w +
+ gradNTT * (-ip_dd.dfT_2.dp_cap_Npart) * Np * w +
// second part of dfT_2/dp_cap
- gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart) * gradNp * w;
+ gradNTT * (-ip_dd.dfT_2.dp_cap_gradNpart) * gradNp * w;
// dfT_2/dT
local_Jac
.template block<temperature_size, temperature_size>(
temperature_index, temperature_index)
- .noalias() -= gradNTT * ip_cv.drho_GR_h_w_eff_dT * NT * w;
+ .noalias() -= gradNTT * ip_dd.dfT_2.dT * NT * w;
// fT_3
- fT.noalias() += NTT *
- (rhoGR * ip_out.darcy_velocity_data.w_GS.transpose() +
- rhoLR * ip_out.darcy_velocity_data.w_LS.transpose()) *
- b * w;
-
- fT.noalias() += gradNTT *
- (current_state.constituent_density_data.rho_C_GR *
- ip_cv.phase_transition_data.hCG *
- ip_out.diffusion_velocity_data.d_CG +
- current_state.constituent_density_data.rho_W_GR *
- ip_cv.phase_transition_data.hWG *
- ip_out.diffusion_velocity_data.d_WG) *
- w;
+ fT.noalias() += NTT * (ip_cv.fT_3.N * w);
+
+ fT.noalias() += gradNTT * ip_cv.fT_3.gradN * w;
// ---------------------------------------------------------------------
// - displacement equation
// ---------------------------------------------------------------------
- KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
+ KUpG.noalias() -= BTI2N * (ip_cv.biot_data() * w);
// dfU_2/dp_GR = dKUpG/dp_GR * p_GR + KUpG. The former is zero, the
// latter is handled below.
- KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w;
+ KUpC.noalias() += BTI2N * (ip_cv.fu_2_KupC.m * w);
// dfU_2/dp_cap = dKUpC/dp_cap * p_cap + KUpC. The former is handled
// here, the latter below.
local_Jac
.template block<displacement_size, W_size>(displacement_index,
W_index)
- .noalias() += BuT * alpha_B * ip_cv.chi_S_L.dchi_dS_L *
- ip_cv.dS_L_dp_cap() * pCap * m * Np * w;
+ .noalias() += BTI2N * (ip_dd.dfu_2_KupC.dp_cap * w);
local_Jac
.template block<displacement_size, displacement_size>(
displacement_index, displacement_index)
- .noalias() += BuT * ip_cd.s_mech_data.stiffness_tensor * Bu * w;
+ .noalias() +=
+ Bu.transpose() * ip_cd.s_mech_data.stiffness_tensor * Bu * w;
// fU_1
- fU.noalias() -= (BuT * current_state.eff_stress_data.sigma -
- N_u_op(Nu).transpose() * rho * b) *
- w;
+ fU.noalias() -=
+ (Bu.transpose() * current_state.eff_stress_data.sigma -
+ N_u_op(Nu).transpose() * ip_cv.volumetric_body_force()) *
+ w;
// KuT
local_Jac
.template block<displacement_size, temperature_size>(
displacement_index, temperature_index)
- .noalias() -=
- BuT *
- (ip_cd.s_mech_data.stiffness_tensor *
- ip_cv.s_therm_exp_data.solid_linear_thermal_expansivity_vector) *
- NT * w;
+ .noalias() -= Bu.transpose() * ip_dd.dfu_1_KuT.dT * NT * w;
/* TODO (naumov) Test with gravity needed to check this Jacobian part.
local_Jac
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index 0d88646ae1fde40a673cc125ef89b63087e69dde..bad9ef32fd5fc973b2f19b53b8b47deff84e9d0c 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -198,6 +198,19 @@ private:
ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const&
models);
+ std::vector<ConstitutiveRelations::DerivativesData<DisplacementDim>>
+ updateConstitutiveVariablesDerivatives(
+ Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev,
+ double const t, double const dt,
+ std::vector<
+ ConstitutiveRelations::ConstitutiveData<DisplacementDim>> const&
+ ip_constitutive_data,
+ std::vector<
+ ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>> const&
+ ip_constitutive_variables,
+ ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const&
+ models);
+
private:
using BMatricesType =
BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>;
diff --git a/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp b/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp
index bc74bda9cf09507481432554aa35bdcc9144697e..7fd5c61d739ecd7064e1400cc7f69407a9691957 100644
--- a/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp
+++ b/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp
@@ -108,7 +108,7 @@ TEST(ProcessLib, TH2MNoPhaseTransition)
rhoWLR);
ASSERT_NEAR(density_water, rhoWLR(), 1e-10);
- CR::EnthalpyData enthalpy;
+ CR::FluidEnthalpyData enthalpy;
CR::MassMoleFractionsData mass_mole_fractions;
CR::FluidDensityData fluid_density;
CR::VapourPartialPressureData vapour_pressure;
diff --git a/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp b/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp
index f636edf49f64ce20807819e0d276115e7285432f..5dcc72fe30c4558ef2e0ac210df31c3ed041271f 100644
--- a/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp
+++ b/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp
@@ -218,7 +218,7 @@ TEST(ProcessLib, TH2MPhaseTransition)
CR::PureLiquidDensityData rhoWLR;
CR::PureLiquidDensityModel rhoWLR_model;
- CR::EnthalpyData enthalpy;
+ CR::FluidEnthalpyData enthalpy;
CR::MassMoleFractionsData mass_mole_fractions;
CR::FluidDensityData fluid_density;
CR::VapourPartialPressureData vapour_pressure;
@@ -480,7 +480,7 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho)
CR::PureLiquidDensityData rhoWLR;
CR::PureLiquidDensityModel rhoWLR_model;
- CR::EnthalpyData enthalpy;
+ CR::FluidEnthalpyData enthalpy;
CR::MassMoleFractionsData mass_mole_fractions;
CR::FluidDensityData fluid_density;
CR::VapourPartialPressureData vapour_pressure;