diff --git a/MaterialLib/MPL/Properties/EffectiveThermalConductivityPorosityMixing.h b/MaterialLib/MPL/Properties/EffectiveThermalConductivityPorosityMixing.h
index f5ef7fdbd78e0d340a293a40dd6bd04bb01f5215..01f80c95962da8ea8bc7e6afe81dead06551a0c8 100644
--- a/MaterialLib/MPL/Properties/EffectiveThermalConductivityPorosityMixing.h
+++ b/MaterialLib/MPL/Properties/EffectiveThermalConductivityPorosityMixing.h
@@ -39,9 +39,9 @@ public:
double const dt) const override;
PropertyDataType dValue(VariableArray const& variable_array,
Variable const variable,
- ParameterLib::SpatialPosition const& /*pos*/,
- double const /*t*/,
- double const /*dt*/) const override;
+ ParameterLib::SpatialPosition const& pos,
+ double const t,
+ double const dt) const override;
private:
ParameterLib::CoordinateSystem const* const local_coordinate_system_;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Base.h b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
index d0f26028eca002d68e9daa55417801abd26c1742..0667e3730d94f05cdd4c2b13063fe78061f6a1e1 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Base.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
@@ -32,6 +32,9 @@ template <int DisplacementDim>
using GlobalDimMatrix =
Eigen::Matrix<double, DisplacementDim, DisplacementDim, Eigen::RowMajor>;
+template <int DisplacementDim>
+using GlobalDimVector = Eigen::Vector<double, DisplacementDim>;
+
struct MediaData
{
explicit MediaData(MaterialPropertyLib::Medium const& medium)
@@ -47,8 +50,8 @@ struct MediaData
struct TemperatureData
{
- double T;
- double T_prev;
+ double T = nan;
+ double T_prev = nan;
};
using ReferenceTemperatureData =
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index 0965f8009c6c145b45c989ca9f5f97f00d7c1fed..7c4611cea705da0a547800e0e7999bc1fd826a7c 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -12,10 +12,13 @@
#include "Biot.h"
#include "Bishops.h"
#include "ConstitutiveDensity.h"
+#include "DarcyVelocity.h"
+#include "DiffusionVelocity.h"
#include "ElasticTangentStiffnessData.h"
#include "Enthalpy.h"
#include "EquivalentPlasticStrainData.h"
#include "FluidDensity.h"
+#include "InternalEnergy.h"
#include "MassMoleFractions.h"
#include "MechanicalStrain.h"
#include "PermeabilityData.h"
@@ -28,9 +31,11 @@
#include "Saturation.h"
#include "SolidCompressibility.h"
#include "SolidDensity.h"
+#include "SolidHeatCapacity.h"
#include "SolidMechanics.h"
#include "SolidThermalExpansion.h"
#include "Swelling.h"
+#include "ThermalConductivity.h"
#include "TotalStress.h"
#include "VapourPartialPressure.h"
#include "Viscosity.h"
@@ -50,6 +55,7 @@ struct StatefulData
MechanicalStrainData<DisplacementDim> mechanical_strain_data;
PureLiquidDensityData rho_W_LR;
ConstituentDensityData constituent_density_data;
+ InternalEnergyData internal_energy_data;
static auto reflect()
{
@@ -70,6 +76,7 @@ struct StatefulDataPrev
PrevState<MechanicalStrainData<DisplacementDim>> mechanical_strain_data;
PrevState<PureLiquidDensityData> rho_W_LR;
PrevState<ConstituentDensityData> constituent_density_data;
+ PrevState<InternalEnergyData> internal_energy_data;
StatefulDataPrev<DisplacementDim>& operator=(
StatefulData<DisplacementDim> const& state)
@@ -80,6 +87,7 @@ struct StatefulDataPrev
mechanical_strain_data = state.mechanical_strain_data;
rho_W_LR = state.rho_W_LR;
constituent_density_data = state.constituent_density_data;
+ internal_energy_data = state.internal_energy_data;
return *this;
}
@@ -97,6 +105,8 @@ struct OutputData
VapourPartialPressureData vapour_pressure_data;
PorosityData porosity_data;
SolidDensityData solid_density_data;
+ DiffusionVelocityData<DisplacementDim> diffusion_velocity_data;
+ DarcyVelocityData<DisplacementDim> darcy_velocity_data;
static auto reflect()
{
@@ -109,7 +119,9 @@ struct OutputData
&Self::fluid_density_data,
&Self::vapour_pressure_data,
&Self::porosity_data,
- &Self::solid_density_data);
+ &Self::solid_density_data,
+ &Self::diffusion_velocity_data,
+ &Self::darcy_velocity_data);
}
};
@@ -138,14 +150,17 @@ struct ConstitutiveTempData
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;
+ EffectiveVolumetricInternalEnergyDerivatives
+ effective_volumetric_internal_energy_d_data;
using DisplacementDimVector = Eigen::Matrix<double, DisplacementDim, 1>;
using DisplacementDimMatrix =
Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
- DisplacementDimMatrix dlambda_dp_GR;
- DisplacementDimMatrix dlambda_dp_cap;
- DisplacementDimMatrix dlambda_dT;
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;
@@ -174,14 +189,6 @@ struct ConstitutiveTempData
DisplacementDimMatrix dadvection_C_dp_cap;
DisplacementDimMatrix dk_over_mu_G_dp_cap;
DisplacementDimMatrix dk_over_mu_L_dp_cap;
- double drho_u_eff_dT = std::numeric_limits<double>::quiet_NaN();
- double drho_u_eff_dp_GR = std::numeric_limits<double>::quiet_NaN();
- double drho_u_eff_dp_cap = std::numeric_limits<double>::quiet_NaN();
- double drho_h_eff_dT = std::numeric_limits<double>::quiet_NaN();
- double drho_h_eff_dp_GR = std::numeric_limits<double>::quiet_NaN();
- double drho_h_eff_dp_cap = std::numeric_limits<double>::quiet_NaN();
- double dh_G_dT = std::numeric_limits<double>::quiet_NaN();
- double dh_L_dT = std::numeric_limits<double>::quiet_NaN();
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();
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index b9526cdf6d34197e8049e328e2930dc1ec84772c..e37d7c906d204a71574007e1dea7db3a946e4bf4 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -20,9 +20,11 @@
#include "Saturation.h"
#include "SolidCompressibility.h"
#include "SolidDensity.h"
+#include "SolidHeatCapacity.h"
#include "SolidMechanics.h"
#include "SolidThermalExpansion.h"
#include "Swelling.h"
+#include "ThermalConductivity.h"
#include "TotalStress.h"
#include "Viscosity.h"
@@ -60,6 +62,7 @@ struct ConstitutiveModels
PermeabilityModel<DisplacementDim> permeability_model;
PureLiquidDensityModel pure_liquid_density_model;
PhaseTransitionModel const& phase_transition_model;
+ ViscosityModel viscosity_model;
#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
PorosityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>
porosity_model;
@@ -69,6 +72,8 @@ struct ConstitutiveModels
PorosityModel porosity_model;
SolidDensityModel solid_density_model;
#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
+ SolidHeatCapacityModel solid_heat_capacity_model;
+ ThermalConductivityModel<DisplacementDim> thermal_conductivity_model;
};
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
new file mode 100644
index 0000000000000000000000000000000000000000..84e66d31d02bcf744f6d5affa6a7183695d2392b
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
@@ -0,0 +1,36 @@
+/**
+ * \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 "ProcessLib/Reflection/ReflectionData.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+struct DarcyVelocityData
+{
+ GlobalDimVector<DisplacementDim> w_GS;
+ GlobalDimVector<DisplacementDim> w_LS;
+
+ static auto reflect()
+ {
+ using Self = DarcyVelocityData<DisplacementDim>;
+ namespace R = ProcessLib::Reflection;
+
+ return std::tuple{
+ R::makeReflectionData("velocity_gas", &Self::w_GS),
+ R::makeReflectionData("velocity_liquid", &Self::w_LS)};
+ }
+};
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h
new file mode 100644
index 0000000000000000000000000000000000000000..25a203ca2c0c03fee94a64527880c3e07f37929e
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.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 "ProcessLib/Reflection/ReflectionData.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+struct DiffusionVelocityData
+{
+ GlobalDimVector<DisplacementDim> d_CG;
+ GlobalDimVector<DisplacementDim> d_WG;
+ GlobalDimVector<DisplacementDim> d_CL;
+ GlobalDimVector<DisplacementDim> d_WL;
+
+ static auto reflect()
+ {
+ using Self = DiffusionVelocityData<DisplacementDim>;
+ namespace R = ProcessLib::Reflection;
+
+ return std::tuple{
+ R::makeReflectionData("diffusion_velocity_gas_gas", &Self::d_CG),
+ R::makeReflectionData("diffusion_velocity_vapour_gas", &Self::d_WG),
+ R::makeReflectionData("diffusion_velocity_solute_liquid",
+ &Self::d_CL),
+ R::makeReflectionData("diffusion_velocity_liquid_liquid",
+ &Self::d_WL)};
+ }
+};
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h
index 48dbbaeb764e5b641ca5f236e9a64e76cea40493..b7a1686618b93faa2386ab2df657d3ae91a92512 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h
@@ -20,6 +20,7 @@ struct EnthalpyData
{
double h_G = nan;
double h_L = nan;
+ double h_S = nan;
static auto reflect()
{
@@ -27,9 +28,22 @@ struct EnthalpyData
namespace R = ProcessLib::Reflection;
return std::tuple{R::makeReflectionData("enthalpy_gas", &Self::h_G),
- R::makeReflectionData("enthalpy_liquid", &Self::h_L)};
+ R::makeReflectionData("enthalpy_liquid", &Self::h_L),
+ R::makeReflectionData("enthalpy_solid", &Self::h_S)};
}
};
+struct EffectiveVolumetricEnthalpy
+{
+ double rho_h_eff = nan;
+};
+
+struct EffectiveVolumetricEnthalpyDerivatives
+{
+ double drho_h_eff_dT = nan;
+ double drho_h_eff_dp_GR = nan;
+ double drho_h_eff_dp_cap = nan;
+};
+
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/FluidDensity.h b/ProcessLib/TH2M/ConstitutiveRelations/FluidDensity.h
index 067bc887273ab6ca921f99dd99cd6b0fd5e34f12..22a0fa60d878c689c00df93cb86e5899d4535845 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/FluidDensity.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/FluidDensity.h
@@ -19,10 +19,10 @@ namespace ConstitutiveRelations
struct FluidDensityData
{
// gas phase density
- double rho_GR = 0.;
+ double rho_GR = nan;
// liquid phase density
- double rho_LR = 0.;
+ double rho_LR = nan;
static auto reflect()
{
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h
new file mode 100644
index 0000000000000000000000000000000000000000..73644146b3452dd53f22f0d4a3b3485d997b0e9a
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h
@@ -0,0 +1,29 @@
+/**
+ * \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 "BaseLib/StrongType.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+using InternalEnergyData =
+ BaseLib::StrongType<double, struct InternalEnergyTag>;
+
+struct EffectiveVolumetricInternalEnergyDerivatives
+{
+ double drho_u_eff_dT = nan;
+ double drho_u_eff_dp_GR = nan;
+ double drho_u_eff_dp_cap = nan;
+};
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/MassMoleFractions.h b/ProcessLib/TH2M/ConstitutiveRelations/MassMoleFractions.h
index a63b2514a47ef459a9e87268578f30b7e0f2b263..a180d179d839ae450dd80709d71a429947eed351 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/MassMoleFractions.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/MassMoleFractions.h
@@ -18,10 +18,10 @@ namespace ConstitutiveRelations
{
struct MassMoleFractionsData
{
- double xnCG = 0.;
- double xmCG = 0.;
- double xnWL = 0.;
- double xmWL = 0.;
+ double xnCG = nan;
+ double xmCG = nan;
+ double xnWL = nan;
+ double xmWL = nan;
static auto reflect()
{
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp
index 9c01e89ce41f1b8991cc05133a25cd688e661cd9..782b8b60ecbec18059c1bd9aef3134bd22c95baa 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp
@@ -43,7 +43,6 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t,
CapillaryPressureData const& p_cap,
TemperatureData const& T_data,
PureLiquidDensityData const& rho_W_LR,
- ViscosityData& viscosity_data,
EnthalpyData& enthalpy_data,
MassMoleFractionsData& mass_mole_fractions_data,
FluidDensityData& fluid_density_data,
@@ -66,10 +65,11 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t,
vapour_pressure_data.pWGR = 0;
// C-component is only component in the gas phase
- cv.xnWG = 0.;
- cv.xmWG = 0.;
- mass_mole_fractions_data.xnCG = 1. - cv.xnWG;
- mass_mole_fractions_data.xmCG = 1. - cv.xmWG;
+ cv.dxmWG_dpGR = 0.;
+ cv.dxmWG_dpCap = 0.;
+ cv.dxmWG_dT = 0.;
+ mass_mole_fractions_data.xnCG = 1.;
+ mass_mole_fractions_data.xmCG = 1.;
auto const M =
gas_phase.property(MaterialPropertyLib::PropertyType::molar_mass)
@@ -80,9 +80,6 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t,
fluid_density_data.rho_GR =
gas_phase.property(MaterialPropertyLib::PropertyType::density)
.template value<double>(variables, x_t.x, x_t.t, x_t.dt);
- viscosity_data.mu_GR =
- gas_phase.property(MaterialPropertyLib::PropertyType::viscosity)
- .template value<double>(variables, x_t.x, x_t.t, x_t.dt);
constituent_density_data.rho_C_GR = fluid_density_data.rho_GR;
constituent_density_data.rho_W_GR = 0;
@@ -90,16 +87,15 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t,
// W-component is only component in the liquid phase
mass_mole_fractions_data.xmWL = 1.;
+ cv.dxmWL_dpGR = 0;
+ cv.dxmWL_dpCap = 0;
+ cv.dxmWL_dT = 0;
auto const pLR = pGR - pCap;
variables.liquid_phase_pressure = pLR;
fluid_density_data.rho_LR = rho_W_LR();
variables.density = fluid_density_data.rho_LR;
- viscosity_data.mu_LR =
- liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity)
- .template value<double>(variables, x_t.x, x_t.t, x_t.dt);
-
// specific heat capacities
auto const cpG =
gas_phase
@@ -121,16 +117,18 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t,
cv.uG = enthalpy_data.h_G - pGR / fluid_density_data.rho_GR;
cv.uL = enthalpy_data.h_L;
- auto const drho_GR_dT =
+ cv.drho_GR_dT =
gas_phase[MaterialPropertyLib::PropertyType::density]
.template dValue<double>(variables,
MaterialPropertyLib::Variable::temperature,
x_t.x, x_t.t, x_t.dt);
- cv.du_G_dT = cpG + pGR * drho_GR_dT / fluid_density_data.rho_GR /
+ cv.du_G_dT = cpG + pGR * cv.drho_GR_dT / fluid_density_data.rho_GR /
fluid_density_data.rho_GR;
cv.du_L_dT = cpL;
+ cv.drho_GR_dp_cap = 0;
+
cv.drho_GR_dp_GR =
gas_phase.property(MaterialPropertyLib::PropertyType::density)
.template dValue<double>(
@@ -141,38 +139,33 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t,
.template dValue<double>(
variables, MaterialPropertyLib::Variable::liquid_phase_pressure,
x_t.x, x_t.t, x_t.dt);
- cv.drho_LR_dp_GR = cv.drho_LR_dp_LR;
cv.du_G_dp_GR = -1 / fluid_density_data.rho_GR +
pGR * cv.drho_GR_dp_GR / fluid_density_data.rho_GR /
fluid_density_data.rho_GR;
cv.drho_C_GR_dp_GR = cv.drho_GR_dp_GR;
- cv.drho_C_LR_dp_LR = 0;
- cv.drho_C_LR_dp_GR = 0;
- cv.drho_C_GR_dT = drho_GR_dT;
+ cv.drho_C_GR_dp_cap = 0;
+ cv.drho_C_GR_dT = cv.drho_GR_dT;
- auto const drho_LR_dT =
+ cv.drho_LR_dT =
liquid_phase[MaterialPropertyLib::PropertyType::density]
.template dValue<double>(variables,
MaterialPropertyLib::Variable::temperature,
x_t.x, x_t.t, x_t.dt);
- cv.drho_C_LR_dT = 0;
-
- cv.du_L_dp_GR = 0;
- cv.du_L_dp_cap = 0;
- /* TODO update to the following when uL has same structure as the uG:
- +-1 / fluid_density_data.rho_LR + pLR* cv.drho_LR_dp_cap /
- fluid_density_data.rho_LR /
- fluid_density_data.rho_LR;
- */
cv.drho_W_LR_dp_LR = cv.drho_LR_dp_LR;
- cv.drho_W_LR_dp_GR = cv.drho_LR_dp_GR;
- cv.drho_W_LR_dT = drho_LR_dT;
+ cv.drho_W_LR_dp_GR = cv.drho_LR_dp_LR;
+ cv.drho_W_LR_dT = cv.drho_LR_dT;
cv.drho_W_GR_dT = 0;
cv.drho_W_GR_dp_GR = 0;
cv.drho_W_GR_dp_cap = 0;
+
+ cv.diffusion_coefficient_vapour = 0.;
+ cv.diffusion_coefficient_solute = 0.;
+
+ cv.hCG = 0;
+ cv.hWG = 0;
}
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h
index 20c7b5abd0ae27dac461618ab3e6551a10c27a0b..918723e906012fb59a5bb1c1cd7c99b67dd23e0c 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,
- ViscosityData& viscosity_data, EnthalpyData& enthalpy_data,
+ EnthalpyData& enthalpy_data,
MassMoleFractionsData& mass_mole_fractions_data,
FluidDensityData& fluid_density_data,
VapourPartialPressureData& vapour_pressure_data,
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PermeabilityData.h b/ProcessLib/TH2M/ConstitutiveRelations/PermeabilityData.h
index b6eeac5c148f8997b6b9c86829efb6f931168679..aea68ad3d26abc9122078c64a0cb2ed2693a33dd 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/PermeabilityData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/PermeabilityData.h
@@ -19,10 +19,10 @@ namespace ConstitutiveRelations
template <int DisplacementDim>
struct PermeabilityData
{
- double k_rel_G;
- double k_rel_L;
- double dk_rel_G_dS_L;
- double dk_rel_L_dS_L;
+ double k_rel_G = nan;
+ double k_rel_L = nan;
+ double dk_rel_G_dS_L = nan;
+ double dk_rel_L_dS_L = nan;
GlobalDimMatrix<DisplacementDim> Ki;
static auto reflect()
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp
index b34faedf64160304bc1dc652a5d4c4221062d314..7f71fcb24920fee9f6fd6094cec574158f19a83d 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp
@@ -112,11 +112,10 @@ PhaseTransition::PhaseTransition(
std::array const required_solvent_component_properties = {
MaterialPropertyLib::specific_heat_capacity};
- std::array const required_gas_properties = {MaterialPropertyLib::density,
- MaterialPropertyLib::viscosity};
+ std::array const required_gas_properties = {MaterialPropertyLib::density};
std::array const required_liquid_properties = {
- MaterialPropertyLib::density, MaterialPropertyLib::viscosity};
+ MaterialPropertyLib::density};
checkRequiredProperties(
gas_phase.component(gas_phase_vapour_component_index_),
@@ -140,7 +139,6 @@ void PhaseTransition::eval(SpaceTimeData const& x_t,
CapillaryPressureData const& p_cap,
TemperatureData const& T_data,
PureLiquidDensityData const& rho_W_LR,
- ViscosityData& viscosity_data,
EnthalpyData& enthalpy_data,
MassMoleFractionsData& mass_mole_fractions_data,
FluidDensityData& fluid_density_data,
@@ -233,9 +231,9 @@ void PhaseTransition::eval(SpaceTimeData const& x_t,
// of the mass balance on the diagonal of the local element
// matrix to be zero). The value is simply made up, seems
// reasonable.
- cv.xnWG =
+ double const xnWG =
std::clamp(vapour_pressure_data.pWGR / pGR, xnWG_min, 1. - xnWG_min);
- mass_mole_fractions_data.xnCG = 1. - cv.xnWG;
+ mass_mole_fractions_data.xnCG = 1. - xnWG;
// gas phase molar fraction derivatives
auto const dxnWG_dpGR = -vapour_pressure_data.pWGR / pGR / pGR;
@@ -243,12 +241,12 @@ void PhaseTransition::eval(SpaceTimeData const& x_t,
auto const dxnWG_dT = dpWGR_dT / pGR;
// molar mass of the gas phase as a mixture of 'air' and vapour
- auto const MG = mass_mole_fractions_data.xnCG * M_C + cv.xnWG * M_W;
+ auto const MG = mass_mole_fractions_data.xnCG * M_C + xnWG * M_W;
variables.molar_mass = MG;
// gas phase mass fractions
- cv.xmWG = cv.xnWG * M_W / MG;
- mass_mole_fractions_data.xmCG = 1. - cv.xmWG;
+ double const xmWG = xnWG * M_W / MG;
+ mass_mole_fractions_data.xmCG = 1. - xmWG;
auto const dxn_dxm_conversion = M_W * M_C / MG / MG;
// gas phase mass fraction derivatives
@@ -302,7 +300,7 @@ void PhaseTransition::eval(SpaceTimeData const& x_t,
// model is still consistent.
constituent_density_data.rho_C_GR =
mass_mole_fractions_data.xmCG * fluid_density_data.rho_GR;
- constituent_density_data.rho_W_GR = cv.xmWG * fluid_density_data.rho_GR;
+ constituent_density_data.rho_W_GR = xmWG * fluid_density_data.rho_GR;
// 'Air'-component partial density derivatives
cv.drho_C_GR_dp_GR = mass_mole_fractions_data.xmCG * cv.drho_GR_dp_GR -
@@ -314,11 +312,11 @@ void PhaseTransition::eval(SpaceTimeData const& x_t,
// Vapour-component partial density derivatives
cv.drho_W_GR_dp_GR =
- cv.xmWG * cv.drho_GR_dp_GR + cv.dxmWG_dpGR * fluid_density_data.rho_GR;
- cv.drho_W_GR_dp_cap = cv.xmWG * cv.drho_GR_dp_cap +
- cv.dxmWG_dpCap * fluid_density_data.rho_GR;
+ xmWG * cv.drho_GR_dp_GR + cv.dxmWG_dpGR * fluid_density_data.rho_GR;
+ cv.drho_W_GR_dp_cap =
+ xmWG * cv.drho_GR_dp_cap + cv.dxmWG_dpCap * fluid_density_data.rho_GR;
cv.drho_W_GR_dT =
- cv.xmWG * cv.drho_GR_dT + cv.dxmWG_dT * fluid_density_data.rho_GR;
+ xmWG * cv.drho_GR_dT + cv.dxmWG_dT * fluid_density_data.rho_GR;
// specific heat capacities of dry air and vapour
auto const cpCG =
@@ -335,11 +333,13 @@ 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 + cv.xmWG * cv.hWG;
+ 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.du_G_dT = 0;
+ cv.du_G_dp_GR = 0;
// diffusion
auto const tortuosity =
@@ -355,11 +355,6 @@ void PhaseTransition::eval(SpaceTimeData const& x_t,
variables.molar_fraction = mass_mole_fractions_data.xnCG;
- // gas phase viscosity
- viscosity_data.mu_GR =
- gas_phase.property(MaterialPropertyLib::PropertyType::viscosity)
- .template value<double>(variables, x_t.x, x_t.t, x_t.dt);
-
// Dissolution part -- Liquid phase properties
// -------------------------------------------
@@ -465,7 +460,7 @@ void PhaseTransition::eval(SpaceTimeData const& x_t,
// liquid phase density derivatives
auto const drhoLR_dpGR = rho_ref_betaP + rho_ref_betaC * dcCL_dpGR;
auto const drhoLR_dpCap = -rho_ref_betaP;
- auto const drhoLR_dT = rho_ref_betaT + rho_ref_betaC * dcCL_dT;
+ cv.drho_LR_dT = rho_ref_betaT + rho_ref_betaC * dcCL_dT;
// solvent partial density derivatives
auto const drhoWLR_dpGR = rho_ref_betaP;
@@ -480,7 +475,7 @@ void PhaseTransition::eval(SpaceTimeData const& x_t,
1. / fluid_density_data.rho_LR *
(drhoWLR_dpCap - mass_mole_fractions_data.xmWL * drhoLR_dpCap);
cv.dxmWL_dT = 1. / fluid_density_data.rho_LR *
- (drhoWLR_dT - mass_mole_fractions_data.xmWL * drhoLR_dT);
+ (drhoWLR_dT - mass_mole_fractions_data.xmWL * cv.drho_LR_dT);
// liquid phase molar fractions and derivatives
mass_mole_fractions_data.xnWL =
@@ -508,12 +503,14 @@ void PhaseTransition::eval(SpaceTimeData const& x_t,
.template value<double>(variables, x_t.x, x_t.t, x_t.dt);
// specific enthalpy of liquid phase and its components
- cv.hCL = cpCL * T + dh_sol;
- cv.hWL = cpWL * T;
- enthalpy_data.h_L = xmCL * cv.hCL + mass_mole_fractions_data.xmWL * cv.hWL;
+ double const hCL = cpCL * T + dh_sol;
+ double const hWL = cpWL * T;
+ 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.du_L_dT = 0;
// diffusion
auto const D_C_L_m =
@@ -522,10 +519,11 @@ void PhaseTransition::eval(SpaceTimeData const& x_t,
cv.diffusion_coefficient_solute =
tortuosity * D_C_L_m; // Note here that D_C_L = D_W_L.
- // liquid phase viscosity
- viscosity_data.mu_LR =
- liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity)
- .template value<double>(variables, x_t.x, x_t.t, x_t.dt);
+ // Some default initializations.
+ cv.drho_LR_dp_LR = 0;
+ cv.drho_W_LR_dp_GR = 0.;
+ cv.drho_W_LR_dT = 0.;
+ cv.drho_W_LR_dp_LR = 0.;
}
} // namespace ConstitutiveRelations
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h
index d4ffbc950f21d2213f4e7acd301782b8bd69b63f..69d621b3dad66159bbf4b350d62e1102a74f893a 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,
- ViscosityData& viscosity_data, EnthalpyData& enthalpy_data,
+ EnthalpyData& 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 2500aa882ae2aead1c12260b741000e4cc878d75..a0042c6afbbad7c47880838830b1ef95681387a5 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h
@@ -17,68 +17,68 @@ namespace ConstitutiveRelations
{
struct PhaseTransitionData
{
- double drho_GR_dp_GR = 0.;
- double drho_GR_dp_cap = 0.;
- double drho_GR_dT = 0.;
+ double drho_GR_dp_GR = nan;
+ double drho_GR_dp_cap = nan;
+ double drho_GR_dT = nan;
- double drho_C_GR_dp_GR = 0.;
- double drho_C_GR_dp_cap = 0.;
- double drho_C_GR_dT = 0.;
+ double drho_C_GR_dp_GR = nan;
+ double drho_C_GR_dp_cap = nan;
+ double drho_C_GR_dT = nan;
- double drho_W_GR_dp_GR = 0.;
- double drho_W_GR_dp_cap = 0.;
- double drho_W_GR_dT = 0.;
+ double drho_W_GR_dp_GR = nan;
+ double drho_W_GR_dp_cap = nan;
+ double drho_W_GR_dT = nan;
- double drho_LR_dp_GR = 0.;
- double drho_LR_dp_cap = 0.;
- double drho_LR_dT = 0.;
- double drho_LR_dp_LR = 0.;
+ double drho_LR_dT = nan;
+ double drho_LR_dp_LR = nan;
- double drho_C_LR_dp_GR = 0.;
- double drho_C_LR_dp_cap = 0.;
- double drho_C_LR_dT = 0.;
- double drho_C_LR_dp_LR = 0.;
+ // TODO (naumov) These three are zero in both models but used in the
+ // assembly. Remove them and simplify assembly or correct the expressions in
+ // the phase transition models?
+ static constexpr double drho_C_LR_dp_GR = 0.;
+ static constexpr double drho_C_LR_dT = 0.;
+ static constexpr double drho_C_LR_dp_LR = 0.;
- double drho_W_LR_dp_GR = 0.;
- double drho_W_LR_dp_cap = 0.;
- double drho_W_LR_dT = 0.;
- double drho_W_LR_dp_LR = 0.;
-
- // constituent mass and molar fractions
- double xnWG = 0.;
- double xmWG = 0.;
+ double drho_W_LR_dp_GR = nan;
+ double drho_W_LR_dT = nan;
+ double drho_W_LR_dp_LR = nan;
// mass fraction derivatives
- double dxmWG_dpGR = 0.;
- double dxmWG_dpCap = 0.;
- double dxmWG_dT = 0.;
+ double dxmWG_dpGR = nan;
+ double dxmWG_dpCap = nan;
+ double dxmWG_dT = nan;
- double dxmWL_dpGR = 0.;
- double dxmWL_dpCap = 0.;
- double dxmWL_dpLR = 0.;
- double dxmWL_dT = 0.;
+ 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 = 0.;
- double diffusion_coefficient_solute = 0.;
+ double diffusion_coefficient_vapour = nan;
+ double diffusion_coefficient_solute = nan;
// specific enthalpies
- double hCG = 0;
- double hWG = 0;
- double hCL = 0;
- double hWL = 0;
+ double hCG = nan;
+ double hWG = nan;
- double dh_G_dT = 0;
- double dh_L_dT = 0;
+ double dh_G_dT = nan;
+ double dh_L_dT = nan;
// specific inner energies
- double uG = 0;
- double uL = 0;
-
- double du_G_dT = 0;
- double du_L_dT = 0;
- double du_G_dp_GR = 0;
- double du_L_dp_GR = 0;
- double du_L_dp_cap = 0;
+ double uG = nan;
+ double uL = nan;
+
+ double du_G_dT = nan;
+ double du_L_dT = nan;
+ double du_G_dp_GR = nan;
+
+ // TODO (naumov) These two are zero in both models but used in the assembly.
+ // Remove them and simplify assembly or correct the expressions in the phase
+ // transition models?
+ static constexpr double du_L_dp_GR = 0;
+ static constexpr double du_L_dp_cap = 0;
};
} // namespace ConstitutiveRelations
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h
index 63672d93664ad04ba7807e5c3ec5dba9f6f2f0ed..8cf92e94f7e03421a82a9c6dd8f7c30178cfe37d 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h
@@ -17,7 +17,6 @@
#include "PhaseTransitionData.h"
#include "PureLiquidDensity.h"
#include "VapourPartialPressure.h"
-#include "Viscosity.h"
namespace ProcessLib::TH2M
{
@@ -33,9 +32,9 @@ struct PhaseTransitionModel
// check for minimum requirement definitions in media object
std::array const required_gas_properties = {
- MaterialPropertyLib::viscosity, MaterialPropertyLib::density};
+ MaterialPropertyLib::density};
std::array const required_liquid_properties = {
- MaterialPropertyLib::viscosity, MaterialPropertyLib::density};
+ MaterialPropertyLib::density};
for (auto const& m : media)
{
@@ -53,7 +52,6 @@ struct PhaseTransitionModel
CapillaryPressureData const& p_cap,
TemperatureData const& T_data,
PureLiquidDensityData const& rho_W_LR,
- ViscosityData& viscosity_data,
EnthalpyData& enthalpy_data,
MassMoleFractionsData& mass_mole_fractions_data,
FluidDensityData& fluid_density_data,
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidHeatCapacity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/SolidHeatCapacity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7e1fcef50cd4f92202ad42ce1d7b93a14b4d08c3
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidHeatCapacity.cpp
@@ -0,0 +1,36 @@
+/**
+ * \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 "SolidHeatCapacity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+
+void SolidHeatCapacityModel::eval(
+ SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ TemperatureData const& T_data,
+ SolidHeatCapacityData& solid_heat_capacity) const
+{
+ namespace MPL = MaterialPropertyLib;
+
+ MPL::VariableArray variables;
+ variables.temperature = T_data.T;
+
+ auto const& mpl_cpS =
+ media_data.solid[MPL::PropertyType::specific_heat_capacity];
+
+ *solid_heat_capacity =
+ mpl_cpS.template value<double>(variables, x_t.x, x_t.t, x_t.dt);
+}
+
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidHeatCapacity.h b/ProcessLib/TH2M/ConstitutiveRelations/SolidHeatCapacity.h
new file mode 100644
index 0000000000000000000000000000000000000000..4f3b4fae0ad1903fac570986ca8e23fa34e746dc
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidHeatCapacity.h
@@ -0,0 +1,32 @@
+/**
+ * \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 "BaseLib/StrongType.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+
+using SolidHeatCapacityData =
+ BaseLib::StrongType<double, struct SolidHeatCapacityTag>;
+
+struct SolidHeatCapacityModel
+{
+ void eval(SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ TemperatureData const& T_data,
+ SolidHeatCapacityData& solid_heat_capacity) const;
+};
+
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h
index 9f225f30f7d2691173ea2bdde69fd1f136db19af..6d4a69a73d1b2d8cfbb5000682858a70cb387d32 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h
@@ -19,8 +19,8 @@ template <int DisplacementDim>
struct SolidThermalExpansionData
{
KelvinVector<DisplacementDim> solid_linear_thermal_expansivity_vector;
- double beta_T_SR; /// Isotropic solid phase volumetric thermal expansion
- /// coefficient.
+ double beta_T_SR = nan; /// Solid phase volumetric thermal expansion
+ /// coefficient.
double thermal_volume_strain = nan;
};
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..abd9e0cc9f3f4d2bdf7cd19f3627209e21b28edd
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp
@@ -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
+ */
+
+#include "ThermalConductivity.h"
+
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+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,
+ ThermalConductivityData<DisplacementDim>& thermal_conductivity_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;
+
+ auto const& mpl_thermal_conductivity =
+ media_data.medium[MPL::PropertyType::thermal_conductivity];
+
+ thermal_conductivity_data.lambda = MPL::formEigenTensor<DisplacementDim>(
+ mpl_thermal_conductivity.value(variables, x_t.x, x_t.t, x_t.dt));
+
+ // Derivatives computed here and not in the MPL property because various
+ // derivatives are not available in the VariableArray.
+
+ auto const lambdaGR =
+ media_data.gas.hasProperty(MPL::PropertyType::thermal_conductivity)
+ ? MPL::formEigenTensor<DisplacementDim>(
+ media_data.gas[MPL::PropertyType::thermal_conductivity].value(
+ variables, x_t.x, x_t.t, x_t.dt))
+ : MPL::formEigenTensor<DisplacementDim>(0.);
+
+ auto const dlambda_GR_dT =
+ media_data.gas.hasProperty(MPL::PropertyType::thermal_conductivity)
+ ? MPL::formEigenTensor<DisplacementDim>(
+ media_data.gas[MPL::PropertyType::thermal_conductivity]
+ .dValue(variables, MPL::Variable::temperature, x_t.x,
+ x_t.t, x_t.dt))
+ : MPL::formEigenTensor<DisplacementDim>(0.);
+
+ auto const lambdaLR =
+ media_data.liquid.hasProperty(MPL::PropertyType::thermal_conductivity)
+ ? MPL::formEigenTensor<DisplacementDim>(
+ media_data.liquid[MPL::PropertyType::thermal_conductivity]
+ .value(variables, x_t.x, x_t.t, x_t.dt))
+ : MPL::formEigenTensor<DisplacementDim>(0.);
+
+ auto const dlambda_LR_dT =
+ media_data.liquid.hasProperty(MPL::PropertyType::thermal_conductivity)
+ ? MPL::formEigenTensor<DisplacementDim>(
+ media_data.liquid[MPL::PropertyType::thermal_conductivity]
+ .dValue(variables, MPL::Variable::temperature, x_t.x,
+ x_t.t, x_t.dt))
+ : MPL::formEigenTensor<DisplacementDim>(0.);
+
+ auto const lambdaSR =
+ media_data.solid.hasProperty(MPL::PropertyType::thermal_conductivity)
+ ? MPL::formEigenTensor<DisplacementDim>(
+ media_data.solid[MPL::PropertyType::thermal_conductivity]
+ .value(variables, x_t.x, x_t.t, x_t.dt))
+ : MPL::formEigenTensor<DisplacementDim>(0.);
+
+ auto const dlambda_SR_dT =
+ media_data.solid.hasProperty(MPL::PropertyType::thermal_conductivity)
+ ? MPL::formEigenTensor<DisplacementDim>(
+ media_data.solid[MPL::PropertyType::thermal_conductivity]
+ .dValue(variables, MPL::Variable::temperature, x_t.x,
+ 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;
+
+ 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_S = 1. - porosity_data.phi;
+
+ // Assuming dS_L/dT = 0, then:
+ // 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 =
+ (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
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h
new file mode 100644
index 0000000000000000000000000000000000000000..fbf34021e55633186144f6fcc267acf6e5d1f45c
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h
@@ -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
+ */
+
+#pragma once
+
+#include "Base.h"
+#include "Porosity.h"
+#include "Saturation.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+
+template <int DisplacementDim>
+struct ThermalConductivityData
+{
+ GlobalDimMatrix<DisplacementDim> lambda;
+ // Currently unused, but there is a comment in TH2MFEM-impl.h referring to
+ // this matrix
+ // GlobalDimMatrix<DisplacementDim> dlambda_dp_GR;
+ GlobalDimMatrix<DisplacementDim> dlambda_dp_cap;
+ GlobalDimMatrix<DisplacementDim> dlambda_dT;
+};
+
+template <int DisplacementDim>
+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;
+};
+
+extern template struct ThermalConductivityModel<2>;
+extern template struct ThermalConductivityModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Viscosity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Viscosity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..782dcb40c5ce7c91f70edcbe476871f441892c41
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Viscosity.cpp
@@ -0,0 +1,38 @@
+/**
+ * \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 "Viscosity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+void ViscosityModel::eval(SpaceTimeData const& x_t, MediaData const& media_data,
+ TemperatureData const& T_data,
+ MassMoleFractionsData const& mass_mole_fractions_data,
+ ViscosityData& viscosity_data) const
+{
+ MaterialPropertyLib::VariableArray variables;
+
+ variables.temperature = T_data.T;
+ variables.molar_fraction = mass_mole_fractions_data.xnCG;
+
+ auto const& liquid_phase = media_data.liquid;
+ auto const& gas_phase = media_data.gas;
+
+ viscosity_data.mu_GR =
+ gas_phase[MaterialPropertyLib::PropertyType::viscosity]
+ .template value<double>(variables, x_t.x, x_t.t, x_t.dt);
+
+ viscosity_data.mu_LR =
+ liquid_phase[MaterialPropertyLib::PropertyType::viscosity]
+ .template value<double>(variables, x_t.x, x_t.t, x_t.dt);
+}
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Viscosity.h b/ProcessLib/TH2M/ConstitutiveRelations/Viscosity.h
index 95a38fd1d24e5aeb68e46916de982a215cae48f2..c46e48c702dc14c6db1e1fe2d9a4c03b8ea41ef3 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Viscosity.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Viscosity.h
@@ -10,6 +10,7 @@
#pragma once
#include "Base.h"
+#include "MassMoleFractions.h"
namespace ProcessLib::TH2M
{
@@ -20,5 +21,13 @@ struct ViscosityData
double mu_GR = nan;
double mu_LR = nan;
};
+
+struct ViscosityModel
+{
+ void eval(SpaceTimeData const& x_t, MediaData const& media_data,
+ TemperatureData const& T_data,
+ MassMoleFractionsData const& mass_mole_fractions_data,
+ ViscosityData& viscosity_data) const;
+};
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h
index ec9de5fdb91b1fc2a3efbb11e14b5b8828cfd52c..0d96e68979e6d82f637ac9f7451d898a56d183ea 100644
--- a/ProcessLib/TH2M/IntegrationPointData.h
+++ b/ProcessLib/TH2M/IntegrationPointData.h
@@ -10,19 +10,12 @@
#pragma once
-#include <memory>
-
-#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
-#include "MathLib/KelvinVector.h"
-#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
-#include "ParameterLib/Parameter.h"
-
namespace ProcessLib
{
namespace TH2M
{
-template <typename BMatricesType, typename ShapeMatrixTypeDisplacement,
- typename ShapeMatricesTypePressure, int DisplacementDim, int NPoints>
+template <typename ShapeMatrixTypeDisplacement,
+ typename ShapeMatricesTypePressure>
struct IntegrationPointData final
{
using GlobalDimMatrixType =
@@ -36,32 +29,7 @@ struct IntegrationPointData final
typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
- // phase enthalpies
- double rho_G_h_G = std::numeric_limits<double>::quiet_NaN();
- double rho_L_h_L = std::numeric_limits<double>::quiet_NaN();
- double rho_S_h_S = std::numeric_limits<double>::quiet_NaN();
-
- // specific enthalpies
- double h_S = std::numeric_limits<double>::quiet_NaN();
-
- // internal energies
- double rho_u_eff = std::numeric_limits<double>::quiet_NaN();
- double rho_u_eff_prev = std::numeric_limits<double>::quiet_NaN();
-
- GlobalDimMatrixType lambda;
- GlobalDimVectorType d_CG;
- GlobalDimVectorType d_WG;
- GlobalDimVectorType d_CL;
- GlobalDimVectorType d_WL;
-
- GlobalDimVectorType w_GS;
- GlobalDimVectorType w_LS;
-
double integration_weight = std::numeric_limits<double>::quiet_NaN();
-
- void pushBackState() { rho_u_eff_prev = rho_u_eff; }
-
- EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
} // namespace TH2M
diff --git a/ProcessLib/TH2M/LocalAssemblerInterface.h b/ProcessLib/TH2M/LocalAssemblerInterface.h
index 9e8d86a0863b988bedb2a1a481c85bd38837a952..30e624d0b542e30c0b3120457bf928a1d196fbb3 100644
--- a/ProcessLib/TH2M/LocalAssemblerInterface.h
+++ b/ProcessLib/TH2M/LocalAssemblerInterface.h
@@ -66,48 +66,6 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
std::string_view name, double const* values,
int const integration_order) = 0;
- virtual std::vector<double> const& getIntPtDarcyVelocityGas(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> const& getIntPtDarcyVelocityLiquid(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> const& getIntPtDiffusionVelocityVapourGas(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> const& getIntPtDiffusionVelocityGasGas(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> const& getIntPtDiffusionVelocitySoluteLiquid(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> const& getIntPtDiffusionVelocityLiquidLiquid(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> const& getIntPtEnthalpySolid(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
// TODO move to NumLib::ExtrapolatableElement
unsigned getNumberOfIntegrationPoints() const
{
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 9203a3d5326e621bb29938ed33525e0f66ee59a6..2f14daaaeebb1d46f078dd6b93b23d245cc5ab5b 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -115,9 +115,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
auto const& medium =
*this->process_data_.media_map.getMedium(this->element_.getID());
- auto const& gas_phase = medium.phase("Gas");
- auto const& liquid_phase = medium.phase("AqueousLiquid");
- auto const& solid_phase = medium.phase("Solid");
ConstitutiveRelations::MediaData media_data{medium};
unsigned const n_integration_points =
@@ -162,7 +159,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ConstitutiveRelations::CapillaryPressureData const pCap_data{pCap};
ConstitutiveRelations::ReferenceTemperatureData const T0{
this->process_data_.reference_temperature(t, pos)[0]};
- double const pLR = pGR - pCap;
GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
GlobalDimVectorType const gradT = gradNp * temperature;
@@ -227,11 +223,15 @@ 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_cv.viscosity_data, ip_out.enthalpy_data,
+ current_state.rho_W_LR, ip_out.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);
+ models.viscosity_model.eval({pos, t, dt}, media_data, T_data,
+ ip_out.mass_mole_fractions_data,
+ ip_cv.viscosity_data);
+
models.porosity_model.eval({pos, t, dt}, media_data,
#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
ip_cv.biot_data, ip_out.eps_data,
@@ -246,20 +246,13 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
ip_out.solid_density_data, ip_cv.solid_density_d_data);
- MPL::VariableArray vars;
- MPL::VariableArray vars_prev;
- vars.temperature = T;
- vars.gas_phase_pressure = pGR;
- vars.capillary_pressure = pCap;
- vars.liquid_phase_pressure = pLR;
-
- // Set volumetric strain for the general case without swelling.
- vars.volumetric_strain = Invariants::trace(eps);
+ models.solid_heat_capacity_model.eval({pos, t, dt}, media_data, T_data,
+ ip_cv.solid_heat_capacity_data);
- vars.liquid_saturation = current_state.S_L_data.S_L;
- vars_prev.liquid_saturation = prev_state.S_L_data->S_L;
-
- vars.porosity = ip_out.porosity_data.phi;
+ 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;
@@ -269,29 +262,20 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
(1. - current_state.S_L_data.S_L) * ip_out.porosity_data.phi;
double const phi_S = 1. - ip_out.porosity_data.phi;
- // thermal conductivity
- ip_data.lambda = MaterialPropertyLib::formEigenTensor<DisplacementDim>(
- medium
- .property(
- MaterialPropertyLib::PropertyType::thermal_conductivity)
- .value(vars, pos, t, dt));
-
- auto const cpS =
- solid_phase.property(MPL::PropertyType::specific_heat_capacity)
- .template value<double>(vars, pos, t, dt);
- ip_data.h_S = cpS * T;
- auto const u_S = ip_data.h_S;
-
- ip_data.rho_u_eff = 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_data.rho_G_h_G =
- phi_G * ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G;
- ip_data.rho_L_h_L =
- phi_L * ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L;
- ip_data.rho_S_h_S =
- phi_S * ip_out.solid_density_data.rho_SR * ip_data.h_S;
+ 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;
@@ -310,45 +294,44 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
// the respective phase transition model, here only the multiplication
// with the gradient of the mass fractions should take place.
- ip_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_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_data.d_WG =
- c.xmWG == 0.
+ 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 / c.xmWG * c.diffusion_coefficient_vapour * gradxmWG;
+ : -phi_G / (1 - ip_out.mass_mole_fractions_data.xmCG) *
+ c.diffusion_coefficient_vapour * gradxmWG;
- ip_data.d_CL =
+ 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_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;
+ 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
// ---------------------------------------------------------------------
- auto const drho_LR_dT =
- liquid_phase.property(MPL::PropertyType::density)
- .template dValue<double>(vars, MPL::Variable::temperature, pos,
- t, dt);
-
- ip_cv.drho_u_eff_dT =
+ 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 * drho_LR_dT * c.uL +
+ 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 * cpS -
+ 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;
@@ -363,58 +346,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
double const dphi_L_dp_cap =
ip_cv.dS_L_dp_cap() * ip_out.porosity_data.phi;
- auto const lambdaGR =
- gas_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
- ? MPL::formEigenTensor<DisplacementDim>(
- gas_phase
- .property(MPL::PropertyType::thermal_conductivity)
- .value(vars, pos, t, dt))
- : MPL::formEigenTensor<DisplacementDim>(0.);
-
- auto const dlambda_GR_dT =
- gas_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
- ? MPL::formEigenTensor<DisplacementDim>(
- gas_phase[MPL::PropertyType::thermal_conductivity].dValue(
- vars, MPL::Variable::temperature, pos, t, dt))
- : MPL::formEigenTensor<DisplacementDim>(0.);
-
- auto const lambdaLR =
- liquid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
- ? MPL::formEigenTensor<DisplacementDim>(
- liquid_phase
- .property(MPL::PropertyType::thermal_conductivity)
- .value(vars, pos, t, dt))
- : MPL::formEigenTensor<DisplacementDim>(0.);
-
- auto const dlambda_LR_dT =
- liquid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
- ? MPL::formEigenTensor<DisplacementDim>(
- liquid_phase[MPL::PropertyType::thermal_conductivity]
- .dValue(vars, MPL::Variable::temperature, pos, t, dt))
- : MPL::formEigenTensor<DisplacementDim>(0.);
-
- auto const lambdaSR =
- solid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
- ? MPL::formEigenTensor<DisplacementDim>(
- solid_phase
- .property(MPL::PropertyType::thermal_conductivity)
- .value(vars, pos, t, dt))
- : MPL::formEigenTensor<DisplacementDim>(0.);
-
- auto const dlambda_SR_dT =
- solid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
- ? MPL::formEigenTensor<DisplacementDim>(
- solid_phase[MPL::PropertyType::thermal_conductivity]
- .dValue(vars, MPL::Variable::temperature, pos, t, dt))
- : MPL::formEigenTensor<DisplacementDim>(0.);
-
- ip_cv.dlambda_dp_cap =
- dphi_G_dp_cap * lambdaGR + dphi_L_dp_cap * lambdaLR;
-
- ip_cv.dlambda_dT = phi_G * dlambda_GR_dT + phi_L * dlambda_LR_dT +
- phi_S * dlambda_SR_dT -
- ip_cv.porosity_d_data.dphi_dT * lambdaSR;
-
// 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)
@@ -423,14 +354,14 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
double const drho_LR_dp_cap = -c.drho_LR_dp_LR;
// drho_GR_dp_cap = 0;
- ip_cv.drho_h_eff_dp_GR =
+ 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.drho_h_eff_dp_cap =
+ 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 +*/
@@ -441,21 +372,23 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
// TODO (naumov) Extend for temperature dependent porosities.
constexpr double dphi_G_dT = 0;
constexpr double dphi_L_dT = 0;
- ip_cv.drho_h_eff_dT =
+ 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 * drho_LR_dT * 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_data.h_S +
- phi_S * ip_cv.solid_density_d_data.drho_SR_dT * ip_data.h_S +
- phi_S * ip_out.solid_density_data.rho_SR * cpS;
+ 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.drho_u_eff_dp_GR =
+ 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 +
@@ -463,7 +396,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
phi_L * drho_LR_dp_GR * c.uL +
phi_L * ip_out.fluid_density_data.rho_LR * c.du_L_dp_GR;
- ip_cv.drho_u_eff_dp_cap =
+ 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 +
@@ -495,14 +428,17 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ip_cv.dS_L_dp_cap() /
ip_cv.viscosity_data.mu_LR;
- ip_data.w_GS = k_over_mu_G * ip_out.fluid_density_data.rho_GR * b -
- k_over_mu_G * gradpGR;
- ip_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_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_data.w_GS +
+ 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 =
@@ -512,7 +448,8 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
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_data.w_LS -
+ -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 =
@@ -521,10 +458,14 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
k_over_mu_L;
ip_cv.drho_GR_h_w_eff_dT =
- c.drho_GR_dT * ip_out.enthalpy_data.h_G * ip_data.w_GS +
- ip_out.fluid_density_data.rho_GR * c.dh_G_dT * ip_data.w_GS +
- drho_LR_dT * ip_out.enthalpy_data.h_L * ip_data.w_LS +
- ip_out.fluid_density_data.rho_LR * c.dh_L_dT * ip_data.w_LS;
+ 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
@@ -968,9 +909,6 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
- auto& ip_data = _ip_data[ip];
- ip_data.pushBackState();
-
this->material_states_[ip].pushBackState();
this->prev_states_[ip] = this->current_states_[ip];
}
@@ -1199,10 +1137,6 @@ void TH2MLocalAssembler<
auto const rho_W_LR_dot =
(current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt;
- auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
-
- auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / 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;
@@ -1388,27 +1322,38 @@ void TH2MLocalAssembler<
// - temperature equation
// ---------------------------------------------------------------------
- MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
+ MTu.noalias() += NTT *
+ ip_cv.effective_volumetric_enthalpy_data.rho_h_eff *
+ mT * Bu * w;
- KTT.noalias() += gradNTT * ip.lambda * gradNT * 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.w_GS +
- rhoLR * ip_out.enthalpy_data.h_L * ip.w_LS) *
+ (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.d_CG +
+ 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.d_WG) *
+ ip_cv.phase_transition_data.hWG *
+ ip_out.diffusion_velocity_data.d_WG) *
w;
- fT.noalias() +=
- NTT * (rhoGR * ip.w_GS.transpose() + rhoLR * ip.w_LS.transpose()) *
- b * w;
+ fT.noalias() += NTT *
+ (rhoGR * ip_out.darcy_velocity_data.w_GS.transpose() +
+ rhoLR * ip_out.darcy_velocity_data.w_LS.transpose()) *
+ b * w;
// ---------------------------------------------------------------------
// - displacement equation
@@ -1690,10 +1635,6 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
auto const rho_W_LR_dot =
(current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt;
- auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
-
- auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / 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;
@@ -2028,30 +1969,39 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
// - temperature equation
// ---------------------------------------------------------------------
- MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
+ MTu.noalias() += NTT *
+ ip_cv.effective_volumetric_enthalpy_data.rho_h_eff *
+ mT * Bu * w;
// dfT_4/dp_GR
// d (MTu * u_dot)/dp_GR
local_Jac
.template block<temperature_size, C_size>(temperature_index,
C_index)
- .noalias() += NTT * ip_cv.drho_h_eff_dp_GR * div_u_dot * NT * w;
+ .noalias() +=
+ NTT * ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_GR *
+ div_u_dot * NT * w;
// dfT_4/dp_cap
// d (MTu * u_dot)/dp_cap
local_Jac
.template block<temperature_size, W_size>(temperature_index,
W_index)
- .noalias() -= NTT * ip_cv.drho_h_eff_dp_cap * div_u_dot * NT * w;
+ .noalias() -=
+ NTT * ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_cap *
+ div_u_dot * NT * w;
// dfT_4/dT
// d (MTu * u_dot)/dT
local_Jac
.template block<temperature_size, temperature_size>(
temperature_index, temperature_index)
- .noalias() += NTT * ip_cv.drho_h_eff_dT * div_u_dot * NT * w;
+ .noalias() +=
+ NTT * ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dT *
+ div_u_dot * NT * w;
- KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
+ KTT.noalias() +=
+ gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
// d KTT/dp_GR * T
// TODO (naumov) always zero if lambda_xR have no derivatives wrt. p_GR.
@@ -2071,40 +2021,57 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
local_Jac
.template block<temperature_size, W_size>(temperature_index,
W_index)
- .noalias() += gradNTT * ip_cv.dlambda_dp_cap * gradT * Np * w;
+ .noalias() += gradNTT *
+ ip_cv.thermal_conductivity_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.dlambda_dT * gradT * NT * w;
+ .noalias() += gradNTT * ip_cv.thermal_conductivity_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;
// dfT_1/dp_GR
local_Jac
.template block<temperature_size, C_size>(temperature_index,
C_index)
- .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_GR * Np * w;
+ .noalias() += NpT * Np *
+ (ip_cv.effective_volumetric_internal_energy_d_data
+ .drho_u_eff_dp_GR /
+ dt * w);
// dfT_1/dp_cap
local_Jac
.template block<temperature_size, W_size>(temperature_index,
W_index)
- .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_cap * Np * w;
+ .noalias() += NpT * Np *
+ (ip_cv.effective_volumetric_internal_energy_d_data
+ .drho_u_eff_dp_cap /
+ dt * w);
// dfT_1/dT
// MTT
local_Jac
.template block<temperature_size, temperature_size>(
temperature_index, temperature_index)
- .noalias() += NTT * ip_cv.drho_u_eff_dT / dt * NT * w;
+ .noalias() +=
+ NTT * NT *
+ (ip_cv.effective_volumetric_internal_energy_d_data.drho_u_eff_dT /
+ dt * w);
// fT_2
fT.noalias() += gradNTT *
- (rhoGR * ip_out.enthalpy_data.h_G * ip.w_GS +
- rhoLR * ip_out.enthalpy_data.h_L * ip.w_LS) *
+ (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;
// dfT_2/dp_GR
@@ -2134,15 +2101,18 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
.noalias() -= gradNTT * ip_cv.drho_GR_h_w_eff_dT * NT * w;
// fT_3
- fT.noalias() +=
- NTT * (rhoGR * ip.w_GS.transpose() + rhoLR * ip.w_LS.transpose()) *
- b * 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() += gradNTT *
(current_state.constituent_density_data.rho_C_GR *
- ip_cv.phase_transition_data.hCG * ip.d_CG +
+ 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.d_WG) *
+ ip_cv.phase_transition_data.hWG *
+ ip_out.diffusion_velocity_data.d_WG) *
w;
// ---------------------------------------------------------------------
@@ -2267,162 +2237,6 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
.noalias() += KUpC;
}
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
- int DisplacementDim>
-std::vector<double> const& TH2MLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
- getIntPtDarcyVelocityGas(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- unsigned const n_integration_points =
- this->integration_method_.getNumberOfPoints();
-
- cache.clear();
- auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
- double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, DisplacementDim, n_integration_points);
-
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- cache_matrix.col(ip).noalias() = _ip_data[ip].w_GS;
- }
-
- return cache;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
- int DisplacementDim>
-std::vector<double> const& TH2MLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
- getIntPtDarcyVelocityLiquid(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- unsigned const n_integration_points =
- this->integration_method_.getNumberOfPoints();
-
- cache.clear();
- auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
- double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, DisplacementDim, n_integration_points);
-
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- cache_matrix.col(ip).noalias() = _ip_data[ip].w_LS;
- }
-
- return cache;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
- int DisplacementDim>
-std::vector<double> const& TH2MLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
- getIntPtDiffusionVelocityVapourGas(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- unsigned const n_integration_points =
- this->integration_method_.getNumberOfPoints();
-
- cache.clear();
- auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
- double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, DisplacementDim, n_integration_points);
-
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- cache_matrix.col(ip).noalias() = _ip_data[ip].d_WG;
- }
-
- return cache;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
- int DisplacementDim>
-std::vector<double> const& TH2MLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
- getIntPtDiffusionVelocityGasGas(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- unsigned const n_integration_points =
- this->integration_method_.getNumberOfPoints();
-
- cache.clear();
- auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
- double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, DisplacementDim, n_integration_points);
-
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- cache_matrix.col(ip).noalias() = _ip_data[ip].d_CG;
- }
-
- return cache;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
- int DisplacementDim>
-std::vector<double> const& TH2MLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
- getIntPtDiffusionVelocitySoluteLiquid(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- unsigned const n_integration_points =
- this->integration_method_.getNumberOfPoints();
-
- cache.clear();
- auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
- double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, DisplacementDim, n_integration_points);
-
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- cache_matrix.col(ip).noalias() = _ip_data[ip].d_CL;
- }
-
- return cache;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
- int DisplacementDim>
-std::vector<double> const& TH2MLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
- getIntPtDiffusionVelocityLiquidLiquid(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- unsigned const n_integration_points =
- this->integration_method_.getNumberOfPoints();
-
- cache.clear();
- auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
- double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, DisplacementDim, n_integration_points);
-
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- cache_matrix.col(ip).noalias() = _ip_data[ip].d_WL;
- }
-
- return cache;
-}
-
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
int DisplacementDim>
void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index da9643ed1706c890f62afe42056776b839226845..0d88646ae1fde40a673cc125ef89b63087e69dde 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -147,7 +147,6 @@ private:
this->solid_material_.initializeInternalStateVariables(
t, x_position, *material_state.material_state_variables);
- ip_data.pushBackState();
material_state.pushBackState();
}
@@ -167,7 +166,6 @@ private:
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
- _ip_data[ip].pushBackState();
this->material_states_[ip].pushBackState();
}
@@ -190,37 +188,6 @@ private:
return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
}
- std::vector<double> const& getIntPtDarcyVelocityGas(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
- std::vector<double> const& getIntPtDarcyVelocityLiquid(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
- std::vector<double> const& getIntPtDiffusionVelocityVapourGas(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
- std::vector<double> const& getIntPtDiffusionVelocityGasGas(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
- std::vector<double> const& getIntPtDiffusionVelocitySoluteLiquid(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
- std::vector<double> const& getIntPtDiffusionVelocityLiquidLiquid(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
-
std::tuple<
std::vector<ConstitutiveRelations::ConstitutiveData<DisplacementDim>>,
std::vector<
@@ -231,28 +198,12 @@ private:
ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const&
models);
- // TODO: Here is some refactoring potential. All secondary variables could
- // be stored in some container to avoid defining one method for each
- // variable.
-
- virtual std::vector<double> const& getIntPtEnthalpySolid(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const override
- {
- return ProcessLib::getIntegrationPointScalarData(_ip_data, &IpData::h_S,
- cache);
- }
-
private:
using BMatricesType =
BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>;
- using IpData =
- IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement,
- ShapeMatricesTypePressure, DisplacementDim,
- ShapeFunctionDisplacement::NPOINTS>;
- std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
+ using IpData = IntegrationPointData<ShapeMatricesTypeDisplacement,
+ ShapeMatricesTypePressure>;
+ std::vector<IpData> _ip_data;
SecondaryData<
typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index 34d2e0d1db4de1e5c312f543bbb073b56b3420b2..b3289ca9d146d002cf4841f474efc27f61ee5d8d 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/ProcessLib/TH2M/TH2MProcess.cpp
@@ -163,33 +163,6 @@ void TH2MProcess<DisplacementDim>::initializeConcreteProcess(
LocalAssemblerInterface<DisplacementDim>::getReflectionDataForOutput(),
_secondary_variables, getExtrapolator(), local_assemblers_);
- add_secondary_variable(
- "velocity_gas", mesh.getDimension(),
- &LocalAssemblerInterface<DisplacementDim>::getIntPtDarcyVelocityGas);
- add_secondary_variable(
- "velocity_liquid", mesh.getDimension(),
- &LocalAssemblerInterface<DisplacementDim>::getIntPtDarcyVelocityLiquid);
- add_secondary_variable(
- "diffusion_velocity_vapour_gas", mesh.getDimension(),
- &LocalAssemblerInterface<
- DisplacementDim>::getIntPtDiffusionVelocityVapourGas);
- add_secondary_variable(
- "diffusion_velocity_gas_gas", mesh.getDimension(),
- &LocalAssemblerInterface<
- DisplacementDim>::getIntPtDiffusionVelocityGasGas);
- add_secondary_variable(
- "diffusion_velocity_solute_liquid", mesh.getDimension(),
- &LocalAssemblerInterface<
- DisplacementDim>::getIntPtDiffusionVelocitySoluteLiquid);
- add_secondary_variable(
- "diffusion_velocity_liquid_liquid", mesh.getDimension(),
- &LocalAssemblerInterface<
- DisplacementDim>::getIntPtDiffusionVelocityLiquidLiquid);
-
- add_secondary_variable(
- "enthalpy_solid", 1,
- &LocalAssemblerInterface<DisplacementDim>::getIntPtEnthalpySolid);
-
//
// enable output of internal variables defined by material models
//
diff --git a/ProcessLib/TH2M/Tests.cmake b/ProcessLib/TH2M/Tests.cmake
index 3c737d92e12425317d86ce084980450f3fef27cf..cb581c572aaff3e86fc4b64ab4233e24bca89972 100644
--- a/ProcessLib/TH2M/Tests.cmake
+++ b/ProcessLib/TH2M/Tests.cmake
@@ -269,7 +269,7 @@ AddTest(
result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu velocity_liquid velocity_liquid 1e-8 1e-8
- result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu sigma sigma 3.3e-6 1e-8
+ result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu sigma sigma 3.4e-6 1e-8
result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu epsilon epsilon 1e-8 1e-8
diff --git a/Tests/Data/TH2M/H2M/OrthotropicSwelling/square.prj b/Tests/Data/TH2M/H2M/OrthotropicSwelling/square.prj
index 83b487a961f428c248d34882ac83b242aca98aec..0726ab51ab76719d871061c94eaf8b4fcfdf136b 100644
--- a/Tests/Data/TH2M/H2M/OrthotropicSwelling/square.prj
+++ b/Tests/Data/TH2M/H2M/OrthotropicSwelling/square.prj
@@ -476,7 +476,7 @@
<vtkdiff>
<regex>square_ts_.*.vtu</regex>
<field>sigma</field>
- <absolute_tolerance>2e-15</absolute_tolerance>
+ <absolute_tolerance>3e-15</absolute_tolerance>
<relative_tolerance>0</relative_tolerance>
</vtkdiff>
<vtkdiff>
diff --git a/Tests/Data/TH2M/HM/Confined_Compression/HM_confined_compression_gas.prj b/Tests/Data/TH2M/HM/Confined_Compression/HM_confined_compression_gas.prj
index 8c95035d44ab31cb0b78b795b5a1743a42f8530c..91676a26672a2a96d461c9c6443099b8ef10a944 100644
--- a/Tests/Data/TH2M/HM/Confined_Compression/HM_confined_compression_gas.prj
+++ b/Tests/Data/TH2M/HM/Confined_Compression/HM_confined_compression_gas.prj
@@ -485,7 +485,7 @@
<vtkdiff>
<regex>HM_confined_compression_gas_ts_.*.vtu</regex>
<field>epsilon</field>
- <absolute_tolerance>5e-15</absolute_tolerance>
+ <absolute_tolerance>6e-15</absolute_tolerance>
<relative_tolerance>1e-15</relative_tolerance>
</vtkdiff>
<vtkdiff>
diff --git a/Tests/Data/TH2M/HM/Confined_Compression/HM_confined_compression_liquid.prj b/Tests/Data/TH2M/HM/Confined_Compression/HM_confined_compression_liquid.prj
index 698a2dfdd75ddd65a9591ff3c739ffc49a80e8c2..9700175157030ca0df9db51ee247a5e72d96de75 100644
--- a/Tests/Data/TH2M/HM/Confined_Compression/HM_confined_compression_liquid.prj
+++ b/Tests/Data/TH2M/HM/Confined_Compression/HM_confined_compression_liquid.prj
@@ -492,7 +492,7 @@
<vtkdiff>
<regex>HM_confined_compression_liquid_ts_.*.vtu</regex>
<field>epsilon</field>
- <absolute_tolerance>5e-15</absolute_tolerance>
+ <absolute_tolerance>6e-15</absolute_tolerance>
<relative_tolerance>1e-15</relative_tolerance>
</vtkdiff>
<vtkdiff>
diff --git a/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp b/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp
index 52b49cd7a029394d5d0aa717258b8c85d238789f..bc74bda9cf09507481432554aa35bdcc9144697e 100644
--- a/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp
+++ b/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp
@@ -34,9 +34,6 @@ TEST(ProcessLib, TH2MNoPhaseTransition)
auto const diffusion_coefficient_vapour = 0.0;
- auto const viscosity_air = 1.e-5;
- auto const viscosity_water = 1.e-3;
-
m << "<medium>\n";
m << " <phases>\n";
@@ -54,7 +51,6 @@ TEST(ProcessLib, TH2MNoPhaseTransition)
// gas phase properties
m << "<properties>\n";
- m << Tests::makeConstantPropertyElement("viscosity", viscosity_air);
m << Tests::makeConstantPropertyElement("density", density_air);
m << Tests::makeConstantPropertyElement("specific_heat_capacity",
specific_heat_capacity_air);
@@ -69,7 +65,6 @@ TEST(ProcessLib, TH2MNoPhaseTransition)
// liquid phase properties
m << "<properties>\n";
- m << Tests::makeConstantPropertyElement("viscosity", viscosity_water);
m << Tests::makeConstantPropertyElement("density", density_water);
m << Tests::makeConstantPropertyElement("specific_heat_capacity",
specific_heat_capacity_water);
@@ -113,7 +108,6 @@ TEST(ProcessLib, TH2MNoPhaseTransition)
rhoWLR);
ASSERT_NEAR(density_water, rhoWLR(), 1e-10);
- CR::ViscosityData viscosity;
CR::EnthalpyData enthalpy;
CR::MassMoleFractionsData mass_mole_fractions;
CR::FluidDensityData fluid_density;
@@ -122,8 +116,8 @@ TEST(ProcessLib, TH2MNoPhaseTransition)
CR::PhaseTransitionData cv;
ptm->eval(x_t, media_data, CR::GasPressureData{pGR},
CR::CapillaryPressureData{pGR}, CR::TemperatureData{T, T}, rhoWLR,
- viscosity, enthalpy, mass_mole_fractions, fluid_density,
- vapour_pressure, constituent_density, cv);
+ enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ constituent_density, cv);
// reference values
double const rhoCGR = density_air;
@@ -141,17 +135,14 @@ TEST(ProcessLib, TH2MNoPhaseTransition)
double const hL = specific_heat_capacity_water * T;
double const uG = hG - pGR / density_air;
double const uL = hL;
- double const muGR = viscosity_air;
- double const muLR = viscosity_water;
ASSERT_NEAR(density_air, fluid_density.rho_GR, 1.0e-10);
ASSERT_NEAR(density_water, fluid_density.rho_LR, 1.0e-10);
ASSERT_NEAR(rhoCGR, constituent_density.rho_C_GR, 1.0e-10);
ASSERT_NEAR(rhoWGR, constituent_density.rho_W_GR, 1.0e-10);
ASSERT_NEAR(rhoCLR, constituent_density.rho_C_LR, 1.0e-10);
- ASSERT_NEAR(xmCG, 1. - cv.xmWG, 1.e-10);
ASSERT_NEAR(xmCG, mass_mole_fractions.xmCG, 1.e-10);
- ASSERT_NEAR(xmWG, cv.xmWG, 1.e-10);
+ ASSERT_NEAR(xmWG, 1 - mass_mole_fractions.xmCG, 1.e-10);
ASSERT_NEAR(dxmWG_dpGR, cv.dxmWG_dpGR, 1.0e-10);
ASSERT_NEAR(dxmCG_dpGR, -cv.dxmWG_dpGR, 1.0e-10);
ASSERT_NEAR(dxmWG_dT, cv.dxmWG_dT, 1.0e-10);
@@ -164,6 +155,4 @@ TEST(ProcessLib, TH2MNoPhaseTransition)
ASSERT_NEAR(uL, cv.uL, 1.0e-10);
ASSERT_NEAR(diffusion_coefficient_vapour, cv.diffusion_coefficient_vapour,
1.0e-10);
- ASSERT_NEAR(muGR, viscosity.mu_GR, 1.0e-10);
- ASSERT_NEAR(muLR, viscosity.mu_LR, 1.0e-10);
}
diff --git a/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp b/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp
index 7fbe098c97de799c3f2068fe39e585a81ff50450..f636edf49f64ce20807819e0d276115e7285432f 100644
--- a/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp
+++ b/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp
@@ -28,8 +28,6 @@ static const double specific_latent_heat_water = 2258000.;
static const double diffusion_coefficient_vapour = 2.6e-6;
static const double thermal_conductivity_air = 0.2;
static const double thermal_conductivity_water = 0.6;
-static const double viscosity_air = 1.e-5;
-static const double viscosity_water = 1.e-3;
// In case of constant gas density
static const double constant_gas_density = 1.;
@@ -100,7 +98,6 @@ std::string MediumDefinition(const bool density_is_constant)
// gas phase properties
m << "<properties>\n";
- m << Tests::makeConstantPropertyElement("viscosity", viscosity_air);
m << Tests::makeConstantPropertyElement("thermal_conductivity",
thermal_conductivity_air);
@@ -149,7 +146,6 @@ std::string MediumDefinition(const bool density_is_constant)
// liquid phase properties
m << "<properties>\n";
- m << Tests::makeConstantPropertyElement("viscosity", viscosity_water);
m << Tests::makeConstantPropertyElement("specific_heat_capacity",
specific_heat_capacity_water);
@@ -222,7 +218,6 @@ TEST(ProcessLib, TH2MPhaseTransition)
CR::PureLiquidDensityData rhoWLR;
CR::PureLiquidDensityModel rhoWLR_model;
- CR::ViscosityData viscosity;
CR::EnthalpyData enthalpy;
CR::MassMoleFractionsData mass_mole_fractions;
CR::FluidDensityData fluid_density;
@@ -246,11 +241,10 @@ TEST(ProcessLib, TH2MPhaseTransition)
rhoWLR_model.eval(x_t, media_data, CR::GasPressureData{pGR + eps_pGR},
p_cap_data, T_data, rhoWLR);
ptm->eval(x_t, media_data, CR::GasPressureData{pGR + eps_pGR},
- p_cap_data, T_data, rhoWLR, viscosity, enthalpy,
- mass_mole_fractions, fluid_density, vapour_pressure,
- constituent_density, cv);
+ p_cap_data, T_data, rhoWLR, enthalpy, mass_mole_fractions,
+ fluid_density, vapour_pressure, constituent_density, cv);
- auto xmWG_plus = cv.xmWG;
+ auto xmWG_plus = 1 - mass_mole_fractions.xmCG;
auto rhoGR_plus = fluid_density.rho_GR;
auto rhoCGR_plus = constituent_density.rho_C_GR;
auto rhoWGR_plus = constituent_density.rho_W_GR;
@@ -258,11 +252,10 @@ TEST(ProcessLib, TH2MPhaseTransition)
rhoWLR_model.eval(x_t, media_data, CR::GasPressureData{pGR - eps_pGR},
p_cap_data, T_data, rhoWLR);
ptm->eval(x_t, media_data, CR::GasPressureData{pGR - eps_pGR},
- p_cap_data, T_data, rhoWLR, viscosity, enthalpy,
- mass_mole_fractions, fluid_density, vapour_pressure,
- constituent_density, cv);
+ p_cap_data, T_data, rhoWLR, enthalpy, mass_mole_fractions,
+ fluid_density, vapour_pressure, constituent_density, cv);
- auto xmWG_minus = cv.xmWG;
+ auto xmWG_minus = 1 - mass_mole_fractions.xmCG;
auto rhoGR_minus = fluid_density.rho_GR;
auto rhoCGR_minus = constituent_density.rho_C_GR;
auto rhoWGR_minus = constituent_density.rho_W_GR;
@@ -271,8 +264,8 @@ TEST(ProcessLib, TH2MPhaseTransition)
rhoWLR_model.eval(x_t, media_data, CR::GasPressureData{pGR}, p_cap_data,
T_data, rhoWLR);
ptm->eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR,
- viscosity, enthalpy, mass_mole_fractions, fluid_density,
- vapour_pressure, constituent_density, cv);
+ enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ constituent_density, cv);
// Central difference derivatives
auto const dxmWG_dpGR = (xmWG_plus - xmWG_minus) / (2. * eps_pGR);
@@ -297,10 +290,10 @@ TEST(ProcessLib, TH2MPhaseTransition)
rhoWLR);
ptm->eval(x_t, media_data, p_GR_data,
CR::CapillaryPressureData{pCap + eps_pCap}, T_data, rhoWLR,
- viscosity, enthalpy, mass_mole_fractions, fluid_density,
- vapour_pressure, constituent_density, cv);
+ enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ constituent_density, cv);
- xmWG_plus = cv.xmWG;
+ xmWG_plus = 1 - mass_mole_fractions.xmCG;
rhoGR_plus = fluid_density.rho_GR;
rhoCGR_plus = constituent_density.rho_C_GR;
rhoWGR_plus = constituent_density.rho_W_GR;
@@ -310,10 +303,10 @@ TEST(ProcessLib, TH2MPhaseTransition)
rhoWLR);
ptm->eval(x_t, media_data, p_GR_data,
CR::CapillaryPressureData{pCap - eps_pCap}, T_data, rhoWLR,
- viscosity, enthalpy, mass_mole_fractions, fluid_density,
- vapour_pressure, constituent_density, cv);
+ enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ constituent_density, cv);
- xmWG_minus = cv.xmWG;
+ xmWG_minus = 1 - mass_mole_fractions.xmCG;
rhoGR_minus = fluid_density.rho_GR;
rhoCGR_minus = constituent_density.rho_C_GR;
rhoWGR_minus = constituent_density.rho_W_GR;
@@ -322,8 +315,8 @@ TEST(ProcessLib, TH2MPhaseTransition)
rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, T_data,
rhoWLR);
ptm->eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR,
- viscosity, enthalpy, mass_mole_fractions, fluid_density,
- vapour_pressure, constituent_density, cv);
+ enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ constituent_density, cv);
// Central difference derivatives
auto const dxmWG_dpCap = (xmWG_plus - xmWG_minus) / (2. * eps_pCap);
@@ -348,11 +341,11 @@ TEST(ProcessLib, TH2MPhaseTransition)
rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data,
CR::TemperatureData{T + eps_T, T}, rhoWLR);
ptm->eval(x_t, media_data, p_GR_data, p_cap_data,
- CR::TemperatureData{T + eps_T, T}, rhoWLR, viscosity,
- enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ CR::TemperatureData{T + eps_T, T}, rhoWLR, enthalpy,
+ mass_mole_fractions, fluid_density, vapour_pressure,
constituent_density, cv);
- xmWG_plus = cv.xmWG;
+ xmWG_plus = 1 - mass_mole_fractions.xmCG;
rhoGR_plus = fluid_density.rho_GR;
rhoCGR_plus = constituent_density.rho_C_GR;
rhoWGR_plus = constituent_density.rho_W_GR;
@@ -360,11 +353,11 @@ TEST(ProcessLib, TH2MPhaseTransition)
rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data,
CR::TemperatureData{T - eps_T, T}, rhoWLR);
ptm->eval(x_t, media_data, p_GR_data, p_cap_data,
- CR::TemperatureData{T - eps_T, T}, rhoWLR, viscosity,
- enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ CR::TemperatureData{T - eps_T, T}, rhoWLR, enthalpy,
+ mass_mole_fractions, fluid_density, vapour_pressure,
constituent_density, cv);
- xmWG_minus = cv.xmWG;
+ xmWG_minus = 1 - mass_mole_fractions.xmCG;
rhoGR_minus = fluid_density.rho_GR;
rhoCGR_minus = constituent_density.rho_C_GR;
rhoWGR_minus = constituent_density.rho_W_GR;
@@ -373,8 +366,8 @@ TEST(ProcessLib, TH2MPhaseTransition)
rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, T_data,
rhoWLR);
ptm->eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR,
- viscosity, enthalpy, mass_mole_fractions, fluid_density,
- vapour_pressure, constituent_density, cv);
+ enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ constituent_density, cv);
// Central difference derivatives
auto const dxmWG_dT = (xmWG_plus - xmWG_minus) / (2. * eps_T);
@@ -397,18 +390,19 @@ TEST(ProcessLib, TH2MPhaseTransition)
// phase pressure
auto const refence_xnWG =
std::clamp(vapour_pressure.pWGR / pGR, 0., 1.);
- ASSERT_NEAR(refence_xnWG, cv.xnWG, 1.e-10);
+ auto const xnWG = 1. - mass_mole_fractions.xnCG;
+ ASSERT_NEAR(refence_xnWG, xnWG, 1.e-10);
// The quotient of constituent partial densities and phase densities
// must be equal to the mass fraction of those constituents in both
// phases.
ASSERT_NEAR(constituent_density.rho_W_GR / fluid_density.rho_GR,
- cv.xmWG, 1.e-10);
+ 1 - mass_mole_fractions.xmCG, 1.e-10);
ASSERT_NEAR(rhoWLR() / fluid_density.rho_LR, mass_mole_fractions.xmWL,
1.e-10);
ASSERT_NEAR(constituent_density.rho_C_GR / fluid_density.rho_GR,
- 1. - cv.xmWG, 1.e-10);
+ mass_mole_fractions.xmCG, 1.e-10);
ASSERT_NEAR(constituent_density.rho_C_LR / fluid_density.rho_LR,
1. - mass_mole_fractions.xmWL, 1.e-10);
@@ -424,8 +418,8 @@ TEST(ProcessLib, TH2MPhaseTransition)
// Gas density (ideal gas in this test):
constexpr double R = MaterialLib::PhysicalConstant::IdealGasConstant;
- auto const xnCG = 1. - cv.xnWG;
- auto const MG = cv.xnWG * molar_mass_water + xnCG * molar_mass_air;
+ auto const MG =
+ xnWG * molar_mass_water + mass_mole_fractions.xnCG * molar_mass_air;
auto const rhoGR = pGR * MG / R / T;
ASSERT_NEAR(rhoGR, fluid_density.rho_GR, 1.e-10);
@@ -449,9 +443,6 @@ TEST(ProcessLib, TH2MPhaseTransition)
diffusion_coefficient_vapour;
ASSERT_NEAR(diffusion_gas_phase, cv.diffusion_coefficient_vapour,
1.0e-10);
-
- ASSERT_NEAR(viscosity_air, viscosity.mu_GR, 1.0e-10);
- ASSERT_NEAR(viscosity_water, viscosity.mu_LR, 1.0e-10);
}
}
@@ -489,7 +480,6 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho)
CR::PureLiquidDensityData rhoWLR;
CR::PureLiquidDensityModel rhoWLR_model;
- CR::ViscosityData viscosity;
CR::EnthalpyData enthalpy;
CR::MassMoleFractionsData mass_mole_fractions;
CR::FluidDensityData fluid_density;
@@ -513,22 +503,20 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho)
rhoWLR_model.eval(x_t, media_data, CR::GasPressureData{pGR + eps_pGR},
p_cap_data, T_data, rhoWLR);
ptm->eval(x_t, media_data, CR::GasPressureData{pGR + eps_pGR},
- p_cap_data, T_data, rhoWLR, viscosity, enthalpy,
- mass_mole_fractions, fluid_density, vapour_pressure,
- constituent_density, cv);
+ p_cap_data, T_data, rhoWLR, enthalpy, mass_mole_fractions,
+ fluid_density, vapour_pressure, constituent_density, cv);
- auto xmWG_plus = cv.xmWG;
+ auto xmWG_plus = 1 - mass_mole_fractions.xmCG;
auto rhoCGR_plus = constituent_density.rho_C_GR;
auto rhoWGR_plus = constituent_density.rho_W_GR;
rhoWLR_model.eval(x_t, media_data, CR::GasPressureData{pGR - eps_pGR},
p_cap_data, T_data, rhoWLR);
ptm->eval(x_t, media_data, CR::GasPressureData{pGR - eps_pGR},
- p_cap_data, T_data, rhoWLR, viscosity, enthalpy,
- mass_mole_fractions, fluid_density, vapour_pressure,
- constituent_density, cv);
+ p_cap_data, T_data, rhoWLR, enthalpy, mass_mole_fractions,
+ fluid_density, vapour_pressure, constituent_density, cv);
- auto xmWG_minus = cv.xmWG;
+ auto xmWG_minus = 1 - mass_mole_fractions.xmCG;
auto rhoCGR_minus = constituent_density.rho_C_GR;
auto rhoWGR_minus = constituent_density.rho_W_GR;
@@ -536,8 +524,8 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho)
rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, T_data,
rhoWLR);
ptm->eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR,
- viscosity, enthalpy, mass_mole_fractions, fluid_density,
- vapour_pressure, constituent_density, cv);
+ enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ constituent_density, cv);
// Central difference derivatives
auto const dxmWG_dpGR = (xmWG_plus - xmWG_minus) / (2. * eps_pGR);
@@ -561,10 +549,10 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho)
rhoWLR);
ptm->eval(x_t, media_data, p_GR_data,
CR::CapillaryPressureData{pCap + eps_pCap}, T_data, rhoWLR,
- viscosity, enthalpy, mass_mole_fractions, fluid_density,
- vapour_pressure, constituent_density, cv);
+ enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ constituent_density, cv);
- xmWG_plus = cv.xmWG;
+ xmWG_plus = 1 - mass_mole_fractions.xmCG;
rhoCGR_plus = constituent_density.rho_C_GR;
rhoWGR_plus = constituent_density.rho_W_GR;
@@ -573,10 +561,10 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho)
rhoWLR);
ptm->eval(x_t, media_data, p_GR_data,
CR::CapillaryPressureData{pCap - eps_pCap}, T_data, rhoWLR,
- viscosity, enthalpy, mass_mole_fractions, fluid_density,
- vapour_pressure, constituent_density, cv);
+ enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ constituent_density, cv);
- xmWG_minus = cv.xmWG;
+ xmWG_minus = 1 - mass_mole_fractions.xmCG;
rhoCGR_minus = constituent_density.rho_C_GR;
rhoWGR_minus = constituent_density.rho_W_GR;
@@ -584,8 +572,8 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho)
rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, T_data,
rhoWLR);
ptm->eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR,
- viscosity, enthalpy, mass_mole_fractions, fluid_density,
- vapour_pressure, constituent_density, cv);
+ enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ constituent_density, cv);
// Central difference derivatives
auto const dxmWG_dpCap = (xmWG_plus - xmWG_minus) / (2. * eps_pCap);
@@ -609,22 +597,22 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho)
rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data,
CR::TemperatureData{T + eps_T, T}, rhoWLR);
ptm->eval(x_t, media_data, p_GR_data, p_cap_data,
- CR::TemperatureData{T + eps_T, T}, rhoWLR, viscosity,
- enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ CR::TemperatureData{T + eps_T, T}, rhoWLR, enthalpy,
+ mass_mole_fractions, fluid_density, vapour_pressure,
constituent_density, cv);
- xmWG_plus = cv.xmWG;
+ xmWG_plus = 1 - mass_mole_fractions.xmCG;
rhoCGR_plus = constituent_density.rho_C_GR;
rhoWGR_plus = constituent_density.rho_W_GR;
rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data,
CR::TemperatureData{T - eps_T, T}, rhoWLR);
ptm->eval(x_t, media_data, p_GR_data, p_cap_data,
- CR::TemperatureData{T - eps_T, T}, rhoWLR, viscosity,
- enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ CR::TemperatureData{T - eps_T, T}, rhoWLR, enthalpy,
+ mass_mole_fractions, fluid_density, vapour_pressure,
constituent_density, cv);
- xmWG_minus = cv.xmWG;
+ xmWG_minus = 1 - mass_mole_fractions.xmCG;
rhoCGR_minus = constituent_density.rho_C_GR;
rhoWGR_minus = constituent_density.rho_W_GR;
@@ -632,8 +620,8 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho)
rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, T_data,
rhoWLR);
ptm->eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR,
- viscosity, enthalpy, mass_mole_fractions, fluid_density,
- vapour_pressure, constituent_density, cv);
+ enthalpy, mass_mole_fractions, fluid_density, vapour_pressure,
+ constituent_density, cv);
// Central difference derivatives
auto const dxmWG_dT = (xmWG_plus - xmWG_minus) / (2. * eps_T);
@@ -655,18 +643,19 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho)
// phase pressure
auto const refence_xnWG =
std::clamp(vapour_pressure.pWGR / pGR, 0., 1.);
- ASSERT_NEAR(refence_xnWG, cv.xnWG, 1.e-10);
+ auto const xnWG = 1. - mass_mole_fractions.xnCG;
+ ASSERT_NEAR(refence_xnWG, xnWG, 1.e-10);
// The quotient of constituent partial densities and phase densities
// must be equal to the mass fraction of those constituents in both
// phases.
ASSERT_NEAR(constituent_density.rho_W_GR / fluid_density.rho_GR,
- cv.xmWG, 1.e-10);
+ 1 - mass_mole_fractions.xmCG, 1.e-10);
ASSERT_NEAR(rhoWLR() / fluid_density.rho_LR, mass_mole_fractions.xmWL,
1.e-10);
ASSERT_NEAR(constituent_density.rho_C_GR / fluid_density.rho_GR,
- 1. - cv.xmWG, 1.e-10);
+ mass_mole_fractions.xmCG, 1.e-10);
ASSERT_NEAR(constituent_density.rho_C_LR / fluid_density.rho_LR,
1. - mass_mole_fractions.xmWL, 1.e-10);
@@ -704,8 +693,5 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho)
diffusion_coefficient_vapour;
ASSERT_NEAR(diffusion_gas_phase, cv.diffusion_coefficient_vapour,
1.0e-10);
-
- ASSERT_NEAR(viscosity_air, viscosity.mu_GR, 1.0e-10);
- ASSERT_NEAR(viscosity_water, viscosity.mu_LR, 1.0e-10);
}
}