diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h
index 3c11736a84ef50fbc2e31a37f0281c3742c8f901..131c83860000b4b68270d841493b8f7feaacda28 100644
--- a/ProcessLib/TH2M/IntegrationPointData.h
+++ b/ProcessLib/TH2M/IntegrationPointData.h
@@ -158,6 +158,8 @@ struct IntegrationPointData final
GlobalDimMatrixType lambda;
GlobalDimVectorType d_CG;
GlobalDimVectorType d_WG;
+ GlobalDimVectorType d_CL;
+ GlobalDimVectorType d_WL;
GlobalDimVectorType w_GS;
GlobalDimVectorType w_LS;
diff --git a/ProcessLib/TH2M/LocalAssemblerInterface.h b/ProcessLib/TH2M/LocalAssemblerInterface.h
index 53a5b52df1de21e3ad24cede912fae49afbd8299..d584d9c734df5dabd16d37d22b1d600f8bbe814a 100644
--- a/ProcessLib/TH2M/LocalAssemblerInterface.h
+++ b/ProcessLib/TH2M/LocalAssemblerInterface.h
@@ -60,6 +60,30 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
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& getIntPtLiquidDensity(
const double t,
std::vector<GlobalVector*> const& x,
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 1cf040c06e4cda4ca73d35d6ee116a6b0e3df788..c2f2ce76d64c572e5db062f5e35e439f50bbd575 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -396,6 +396,7 @@ std::vector<ConstitutiveVariables<DisplacementDim>> TH2MLocalAssembler<
ip_data.xmCG = 1. - c.xmWG;
ip_data.xmWG = c.xmWG;
ip_data.xmWL = c.xmWL;
+ auto const xmCL = 1. - c.xmWL;
ip_data.diffusion_coefficient_vapour = c.diffusion_coefficient_vapour;
ip_data.diffusion_coefficient_solute = c.diffusion_coefficient_solute;
@@ -411,6 +412,11 @@ std::vector<ConstitutiveVariables<DisplacementDim>> TH2MLocalAssembler<
ip_data.dxmWG_dT * gradT;
const GlobalDimVectorType gradxmCG = -gradxmWG;
+ const GlobalDimVectorType gradxmWL = ip_data.dxmWL_dpGR * gradpGR +
+ ip_data.dxmWL_dpCap * gradpCap +
+ ip_data.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.
@@ -429,6 +435,19 @@ std::vector<ConstitutiveVariables<DisplacementDim>> TH2MLocalAssembler<
ip_data.diffusion_coefficient_vapour *
gradxmWG;
+ ip_data.d_CL = xmCL == 0. ? 0. * gradxmCL // Keep d_CL's dimension and
+ // prevent division by zero
+ : -phi_L / xmCL *
+ ip_data.diffusion_coefficient_solute *
+ gradxmCL;
+
+ ip_data.d_WL = ip_data.xmWL == 0.
+ ? 0. * gradxmWL // Keep d_WG's dimension and prevent
+ // division by zero
+ : -phi_L / ip_data.xmWL *
+ ip_data.diffusion_coefficient_solute *
+ gradxmWL;
+
// ---------------------------------------------------------------------
// Derivatives for Jacobian
// ---------------------------------------------------------------------
@@ -2154,6 +2173,110 @@ std::vector<double> const& TH2MLocalAssembler<
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 =
+ _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 =
+ _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 =
+ _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 =
+ _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 8691f2e48fc8a0eb976c0490f89c0a185aebd310..4535cd2f03d45c363773fae001f6cfff97242235 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -179,6 +179,26 @@ private:
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::vector<ConstitutiveVariables<DisplacementDim>>
updateConstitutiveVariables(Eigen::VectorXd const& local_x,
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index bf2a500474671ea9f645865910d955c61bfda901..2ba3ddf98d6f7d1bbc5ee8f471ca15c099ed84a8 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/ProcessLib/TH2M/TH2MProcess.cpp
@@ -203,6 +203,19 @@ void TH2MProcess<DisplacementDim>::initializeConcreteProcess(
add_secondary_variable(
"velocity_liquid", mesh.getDimension(),
&LocalAssemblerInterface::getIntPtDarcyVelocityLiquid);
+ add_secondary_variable(
+ "diffusion_velocity_vapour_gas", mesh.getDimension(),
+ &LocalAssemblerInterface::getIntPtDiffusionVelocityVapourGas);
+ add_secondary_variable(
+ "diffusion_velocity_gas_gas", mesh.getDimension(),
+ &LocalAssemblerInterface::getIntPtDiffusionVelocityGasGas);
+ add_secondary_variable(
+ "diffusion_velocity_solute_liquid", mesh.getDimension(),
+ &LocalAssemblerInterface::getIntPtDiffusionVelocitySoluteLiquid);
+ add_secondary_variable(
+ "diffusion_velocity_liquid_liquid", mesh.getDimension(),
+ &LocalAssemblerInterface::getIntPtDiffusionVelocityLiquidLiquid);
+
add_secondary_variable("saturation", 1,
&LocalAssemblerInterface::getIntPtSaturation);
add_secondary_variable("vapour_pressure", 1,