diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Base.h b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
index e5ec55f00b49947babc5bbb64c95ea1a5de9dd3b..70c0467e6cf95b44b5e6bd4da68ca8f3b574574e 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Base.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
@@ -67,6 +67,10 @@ template <int DisplacementDim>
using CapillaryPressureGradientData =
BaseLib::StrongType<GlobalDimVector<DisplacementDim>,
struct CapillaryPressureGradientTag>;
+template <int DisplacementDim>
+using TemperatureGradientData =
+ BaseLib::StrongType<GlobalDimVector<DisplacementDim>,
+ struct TemperatureGradientTag>;
template <int DisplacementDim>
using SpecificBodyForceData =
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index 3ac2a4840765c8fa7bf730eb159cc3c0a068ed78..6c80e242c27b7ef77cf90d27e8b3c1e5dbb34d53 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -13,6 +13,7 @@
#include "Biot.h"
#include "Bishops.h"
#include "CEquation.h"
+#include "DiffusionVelocity.h"
#include "ElasticTangentStiffnessModel.h"
#include "Enthalpy.h"
#include "Gravity.h"
@@ -85,6 +86,7 @@ struct ConstitutiveModels
InternalEnergyModel internal_energy_model;
EffectiveVolumetricEnthalpyModel effective_volumetric_enthalpy_model;
GravityModel<DisplacementDim> gravity_model;
+ DiffusionVelocityModel<DisplacementDim> diffusion_velocity_model;
DarcyVelocityModel<DisplacementDim> darcy_velocity_model;
FC1Model<DisplacementDim> fC_1_model;
FC2aModel fC_2a_model;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..20bd591702698b46790543ed7d87a73b4512a6d0
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.cpp
@@ -0,0 +1,90 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "DiffusionVelocity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void DiffusionVelocityModel<DisplacementDim>::eval(
+ CapillaryPressureGradientData<DisplacementDim> const& grad_p_cap,
+ GasPressureGradientData<DisplacementDim> const& grad_p_GR,
+ MassMoleFractionsData const& mass_mole_fractions_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ TemperatureGradientData<DisplacementDim> const& grad_T,
+ DiffusionVelocityData<DisplacementDim>& diffusion_velocity_data) const
+{
+ auto const gradxmWG = phase_transition_data.dxmWG_dpGR * grad_p_GR() +
+ phase_transition_data.dxmWG_dpCap * grad_p_cap() +
+ phase_transition_data.dxmWG_dT * grad_T();
+ auto const gradxmCG = -gradxmWG;
+
+ auto const gradxmWL = phase_transition_data.dxmWL_dpGR * grad_p_GR() +
+ phase_transition_data.dxmWL_dpCap * grad_p_cap() +
+ phase_transition_data.dxmWL_dT * grad_T();
+ auto const gradxmCL = -gradxmWL;
+
+ double const phi_L = S_L_data.S_L * porosity_data.phi;
+ double const phi_G = (1. - S_L_data.S_L) * porosity_data.phi;
+
+ double const phi_G_D_vapour =
+ phi_G * phase_transition_data.diffusion_coefficient_vapour;
+ double const phi_L_D_solute =
+ phi_L * phase_transition_data.diffusion_coefficient_solute;
+
+ if (mass_mole_fractions_data.xmCG == 0)
+ {
+ diffusion_velocity_data.d_CG = GlobalDimVector<DisplacementDim>::Zero();
+ }
+ else
+ {
+ diffusion_velocity_data.d_CG =
+ -phi_G_D_vapour / mass_mole_fractions_data.xmCG * gradxmCG;
+ }
+
+ if (mass_mole_fractions_data.xmCG == 1)
+ {
+ diffusion_velocity_data.d_WG = GlobalDimVector<DisplacementDim>::Zero();
+ }
+ else
+ {
+ diffusion_velocity_data.d_WG =
+ -phi_G_D_vapour / (1 - mass_mole_fractions_data.xmCG) * gradxmWG;
+ }
+
+ if (mass_mole_fractions_data.xmWL == 1)
+
+ {
+ diffusion_velocity_data.d_CL = GlobalDimVector<DisplacementDim>::Zero();
+ }
+ else
+ {
+ diffusion_velocity_data.d_CL =
+ -phi_L_D_solute / (1. - mass_mole_fractions_data.xmWL) * gradxmCL;
+ }
+
+ if (mass_mole_fractions_data.xmWL == 0)
+ {
+ diffusion_velocity_data.d_WL = GlobalDimVector<DisplacementDim>::Zero();
+ }
+ else
+ {
+ diffusion_velocity_data.d_WL =
+ -phi_L_D_solute / mass_mole_fractions_data.xmWL * gradxmWL;
+ }
+}
+
+template struct DiffusionVelocityModel<2>;
+template struct DiffusionVelocityModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h
index 25a203ca2c0c03fee94a64527880c3e07f37929e..7ec45405ed82137de5b2d6ae8206c1ad10ff5833 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h
@@ -10,7 +10,11 @@
#pragma once
#include "Base.h"
+#include "MassMoleFractions.h"
+#include "PhaseTransitionData.h"
+#include "Porosity.h"
#include "ProcessLib/Reflection/ReflectionData.h"
+#include "Saturation.h"
namespace ProcessLib::TH2M
{
@@ -38,5 +42,22 @@ struct DiffusionVelocityData
&Self::d_WL)};
}
};
+
+template <int DisplacementDim>
+struct DiffusionVelocityModel
+{
+ void eval(
+ CapillaryPressureGradientData<DisplacementDim> const& grad_p_cap,
+ GasPressureGradientData<DisplacementDim> const& grad_p_GR,
+ MassMoleFractionsData const& mass_mole_fractions_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ TemperatureGradientData<DisplacementDim> const& grad_T,
+ DiffusionVelocityData<DisplacementDim>& diffusion_velocity_data) const;
+};
+
+extern template struct DiffusionVelocityModel<2>;
+extern template struct DiffusionVelocityModel<3>;
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index c87165fd500875929b65544adb9737859ba26bd7..f354e92facdae518d459ea69330aca4f83e3893f 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -166,6 +166,8 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
grad_p_GR{gradpGR};
ConstitutiveRelations::CapillaryPressureGradientData<
DisplacementDim> const grad_p_cap{gradpCap};
+ ConstitutiveRelations::TemperatureGradientData<DisplacementDim> const
+ grad_T{gradT};
// medium properties
models.elastic_tangent_stiffness_model.eval({pos, t, dt}, T_data,
@@ -272,12 +274,14 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
this->process_data_.specific_body_force},
ip_cv.volumetric_body_force);
- auto const& c = ip_cv.phase_transition_data;
-
- auto const phi_L =
- current_state.S_L_data.S_L * ip_out.porosity_data.phi;
- auto const phi_G =
- (1. - current_state.S_L_data.S_L) * ip_out.porosity_data.phi;
+ models.diffusion_velocity_model.eval(grad_p_cap,
+ grad_p_GR,
+ ip_out.mass_mole_fractions_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ grad_T,
+ ip_out.diffusion_velocity_data);
ip_out.enthalpy_data.h_S = ip_cv.solid_heat_capacity_data() * T;
@@ -443,50 +447,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
prev_state.internal_energy_data,
ip_cv.fT_1);
- // for variable output
- auto const xmCL = 1. - ip_out.mass_mole_fractions_data.xmWL;
-
- const GlobalDimVectorType gradxmWG = c.dxmWG_dpGR * gradpGR +
- c.dxmWG_dpCap * gradpCap +
- c.dxmWG_dT * gradT;
- const GlobalDimVectorType gradxmCG = -gradxmWG;
-
- const GlobalDimVectorType gradxmWL = c.dxmWL_dpGR * gradpGR +
- c.dxmWL_dpCap * gradpCap +
- c.dxmWL_dT * gradT;
- const GlobalDimVectorType gradxmCL = -gradxmWL;
-
- // Todo: factor -phiAlpha / xmZetaAlpha * DZetaAlpha can be evaluated in
- // the respective phase transition model, here only the multiplication
- // with the gradient of the mass fractions should take place.
-
- ip_out.diffusion_velocity_data.d_CG =
- ip_out.mass_mole_fractions_data.xmCG == 0.
- ? 0. * gradxmCG // Keep d_CG's dimension and prevent
- // division by zero
- : -phi_G / ip_out.mass_mole_fractions_data.xmCG *
- c.diffusion_coefficient_vapour * gradxmCG;
-
- ip_out.diffusion_velocity_data.d_WG =
- ip_out.mass_mole_fractions_data.xmCG == 1.
- ? 0. * gradxmWG // Keep d_WG's dimension and prevent
- // division by zero
- : -phi_G / (1 - ip_out.mass_mole_fractions_data.xmCG) *
- c.diffusion_coefficient_vapour * gradxmWG;
-
- ip_out.diffusion_velocity_data.d_CL =
- xmCL == 0.
- ? 0. * gradxmCL // Keep d_CL's dimension and
- // prevent division by zero
- : -phi_L / xmCL * c.diffusion_coefficient_solute * gradxmCL;
-
- ip_out.diffusion_velocity_data.d_WL =
- ip_out.mass_mole_fractions_data.xmWL == 0.
- ? 0. * gradxmWL // Keep d_WG's dimension and prevent
- // division by zero
- : -phi_L / ip_out.mass_mole_fractions_data.xmWL *
- c.diffusion_coefficient_solute * gradxmWL;
-
// ---------------------------------------------------------------------
// Derivatives for Jacobian
// ---------------------------------------------------------------------