diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index fc9b04d0feab300e5f9db3e4edb6fee908d3a971..3860d4374c6ba873623ea1f3fd881690681b48b2 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -171,6 +171,7 @@ struct ConstitutiveTempData
FW2Data fW_2;
FW3aData fW_3a;
FW4LWpGData<DisplacementDim> fW_4_LWpG;
+ FW4LWpCData<DisplacementDim> fW_4_LWpC;
using DisplacementDimVector = Eigen::Matrix<double, DisplacementDim, 1>;
using DisplacementDimMatrix =
@@ -181,13 +182,6 @@ struct ConstitutiveTempData
DisplacementDimVector drho_LR_h_w_eff_dp_cap_Npart;
DisplacementDimMatrix drho_LR_h_w_eff_dp_cap_gradNpart;
DisplacementDimVector drho_GR_h_w_eff_dT;
- DisplacementDimMatrix dfW_4_LWpC_a_dp_GR;
- DisplacementDimMatrix dfW_4_LWpC_a_dp_cap;
- DisplacementDimMatrix dfW_4_LWpC_a_dT;
- DisplacementDimMatrix dfW_4_LWpC_d_dp_GR;
- DisplacementDimMatrix dfW_4_LWpC_d_dp_cap;
- DisplacementDimMatrix dfW_4_LWpC_d_dT;
- DisplacementDimMatrix dk_over_mu_L_dp_cap;
};
/// Data that stores intermediate values (derivatives), which are not needed
@@ -212,6 +206,7 @@ struct DerivativesData
FW2DerivativeData dfW_2;
FW3aDerivativeData dfW_3a;
FW4LWpGDerivativeData<DisplacementDim> dfW_4_LWpG;
+ FW4LWpCDerivativeData<DisplacementDim> dfW_4_LWpC;
};
} // namespace ConstitutiveRelations
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index 33c1eaf495e32b3c3ba0b128d0fe2a0cf2990ab4..3f0ee616b7b4b9a0a8b0ec22d57183d482e2c920 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -97,6 +97,7 @@ struct ConstitutiveModels
FW2Model fW_2_model;
FW3aModel fW_3a_model;
FW4LWpGModel<DisplacementDim> fW_4_LWpG_model;
+ FW4LWpCModel<DisplacementDim> fW_4_LWpC_model;
};
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
index dc127ce8621f132d40464d0b728d17106aed6914..27495f7dbe0892ecd2ceb49f15272753b747c179 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
@@ -254,5 +254,96 @@ void FW4LWpGModel<DisplacementDim>::dEval(
template struct FW4LWpGModel<2>;
template struct FW4LWpGModel<3>;
+template <int DisplacementDim>
+void FW4LWpCModel<DisplacementDim>::eval(
+ AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FW4LWpCData<DisplacementDim>& fW_4_LWpC) const
+{
+ double const sD_G = phase_transition_data.diffusion_coefficient_vapour;
+ double const sD_L = phase_transition_data.diffusion_coefficient_solute;
+
+ double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi;
+ double const phi_L = S_L_data.S_L * porosity_data.phi;
+
+ double const diffusion_WGpCap = phi_G * fluid_density_data.rho_GR * sD_G *
+ phase_transition_data.dxmWG_dpCap;
+ double const diffusion_WLpCap = phi_L * fluid_density_data.rho_LR * sD_L *
+ phase_transition_data.dxmWL_dpCap;
+
+ double const diffusion_W_pCap = diffusion_WGpCap + diffusion_WLpCap;
+
+ auto const I =
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
+
+ fW_4_LWpC.L.noalias() = diffusion_W_pCap * I - advection_data.advection_W_L;
+
+ // Jac wrong?
+ // double const diffusion_W_L_p = phi_L * fluid_density_data.rho_LR * sD_L *
+ // phase_transition_data.dxmWL_dpLR;
+ //
+ // fW_4_LWpC.L.noalias() = diffusion_W_L_p * I +
+ // advection_data.advection_W_L;
+}
+
+template <int DisplacementDim>
+void FW4LWpCModel<DisplacementDim>::dEval(
+ AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ ViscosityData const& viscosity_data,
+ FW4LWpCDerivativeData<DisplacementDim>& dfW_4_LWpC) const
+{
+ ////// Diffusion Part /////
+ // TODO (naumov) d(diffusion_W_pCap)/dX for dxmW*/d* != 0
+
+ ////// Advection part /////
+ GlobalDimMatrix<DisplacementDim> const k_over_mu_L =
+ permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR;
+
+ dfW_4_LWpC.dp_GR = phase_transition_data.drho_W_LR_dp_GR * k_over_mu_L
+ //+ rhoWLR * (dk_over_mu_L_dp_GR = 0)
+ ;
+
+ double const sD_G = phase_transition_data.diffusion_coefficient_vapour;
+ double const sD_L = phase_transition_data.diffusion_coefficient_solute;
+
+ double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi;
+ double const phi_L = S_L_data.S_L * porosity_data.phi;
+
+ double const diffusion_WGpCap = phi_G * fluid_density_data.rho_GR * sD_G *
+ phase_transition_data.dxmWG_dpCap;
+ double const diffusion_WLpCap = phi_L * fluid_density_data.rho_LR * sD_L *
+ phase_transition_data.dxmWL_dpCap;
+
+ double const diffusion_W_pCap = diffusion_WGpCap + diffusion_WLpCap;
+
+ auto const I =
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
+
+ dfW_4_LWpC.dp_cap = diffusion_W_pCap * I - advection_data.advection_W_L;
+
+ auto const dk_over_mu_L_dp_cap = permeability_data.Ki *
+ permeability_data.dk_rel_L_dS_L *
+ dS_L_dp_cap() / viscosity_data.mu_LR;
+ dfW_4_LWpC.dp_cap = -phase_transition_data.drho_W_LR_dp_LR * k_over_mu_L +
+ rho_W_LR() * dk_over_mu_L_dp_cap;
+
+ dfW_4_LWpC.dT = phase_transition_data.drho_W_LR_dT * k_over_mu_L
+ //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T))
+ ;
+}
+
+template struct FW4LWpCModel<2>;
+template struct FW4LWpCModel<3>;
+
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
index 8d88c9e7eff26175c195dc1171de48325482ba7a..206d3fedafe237263956d7b3405a534742a82505 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
@@ -149,5 +149,44 @@ struct FW4LWpGModel
extern template struct FW4LWpGModel<2>;
extern template struct FW4LWpGModel<3>;
+template <int DisplacementDim>
+struct FW4LWpCData
+{
+ GlobalDimMatrix<DisplacementDim> L;
+};
+
+template <int DisplacementDim>
+struct FW4LWpCDerivativeData
+{
+ GlobalDimMatrix<DisplacementDim> dp_GR;
+ GlobalDimMatrix<DisplacementDim> dp_cap;
+ GlobalDimMatrix<DisplacementDim> dT;
+};
+
+template <int DisplacementDim>
+struct FW4LWpCModel
+{
+ void eval(AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ FW4LWpCData<DisplacementDim>& fW_4_LWpC) const;
+
+ void dEval(AdvectionData<DisplacementDim> const& advection_data,
+ FluidDensityData const& fluid_density_data,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ PhaseTransitionData const& phase_transition_data,
+ PorosityData const& porosity_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ ViscosityData const& viscosity_data,
+ FW4LWpCDerivativeData<DisplacementDim>& dfW_4_LWpC) const;
+};
+
+extern template struct FW4LWpCModel<2>;
+extern template struct FW4LWpCModel<3>;
+
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 934bf873412441bcc33dfc227d16452b63a05c26..394beeabb9b14b173df127d51426cb4b67bb9c59 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -380,6 +380,13 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
current_state.S_L_data,
ip_cv.fW_4_LWpG);
+ models.fW_4_LWpC_model.eval(ip_cv.advection_data,
+ ip_out.fluid_density_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_cv.fW_4_LWpC);
+
// for variable output
auto const xmCL = 1. - ip_out.mass_mole_fractions_data.xmWL;
@@ -436,11 +443,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
ip_cv.viscosity_data.mu_LR;
- ip_cv.dk_over_mu_L_dp_cap = ip_out.permeability_data.Ki *
- ip_out.permeability_data.dk_rel_L_dS_L *
- ip_cv.dS_L_dp_cap() /
- ip_cv.viscosity_data.mu_LR;
-
ip_out.darcy_velocity_data.w_GS =
k_over_mu_G * ip_out.fluid_density_data.rho_GR * b -
k_over_mu_G * gradpGR;
@@ -487,24 +489,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
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
-
- ip_cv.dfW_4_LWpC_a_dp_GR = c.drho_W_LR_dp_GR * k_over_mu_L
- //+ rhoWLR * (dk_over_mu_L_dp_GR = 0)
- ;
- ip_cv.dfW_4_LWpC_a_dp_cap =
- -c.drho_W_LR_dp_LR * k_over_mu_L +
- current_state.rho_W_LR() * ip_cv.dk_over_mu_L_dp_cap;
- ip_cv.dfW_4_LWpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
- //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T))
- ;
-
- // TODO (naumov) for dxmW*/d* != 0
- ip_cv.dfW_4_LWpC_d_dp_GR =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
- ip_cv.dfW_4_LWpC_d_dp_cap =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
- ip_cv.dfW_4_LWpC_d_dT =
- Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
}
return {ip_constitutive_data, ip_constitutive_variables};
@@ -708,6 +692,17 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ip_cv.dS_L_dp_cap,
ip_cv.viscosity_data,
ip_dd.dfW_4_LWpG);
+
+ models.fW_4_LWpC_model.dEval(ip_cv.advection_data,
+ ip_out.fluid_density_data,
+ ip_out.permeability_data,
+ ip_cv.phase_transition_data,
+ ip_out.porosity_data,
+ current_state.rho_W_LR,
+ current_state.S_L_data,
+ ip_cv.dS_L_dp_cap,
+ ip_cv.viscosity_data,
+ ip_dd.dfW_4_LWpC);
}
return ip_d_data;
@@ -1177,27 +1172,17 @@ void TH2MLocalAssembler<
using DisplacementDimMatrix =
Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
- DisplacementDimMatrix const diffusion_WGpCap =
- phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpCap;
- DisplacementDimMatrix const diffusion_WLpCap =
- phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpCap;
-
DisplacementDimMatrix const diffusion_WGT =
phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dT;
DisplacementDimMatrix const diffusion_WLT =
phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dT;
- DisplacementDimMatrix const diffusion_W_pCap =
- diffusion_WGpCap + diffusion_WLpCap;
-
DisplacementDimMatrix const diffusion_W_T =
diffusion_WGT + diffusion_WLT;
LWpG.noalias() += gradNpT * ip_cv.fW_4_LWpG.L * gradNp * w;
- LWpC.noalias() +=
- gradNpT * (diffusion_W_pCap - ip_cv.advection_data.advection_W_L) *
- gradNp * w;
+ LWpC.noalias() += gradNpT * ip_cv.fW_4_LWpC.L * gradNp * w;
LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
@@ -1665,8 +1650,6 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
- GlobalDimMatrixType const diffusion_W_L_p =
- phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpLR;
GlobalDimMatrixType const diffusion_W_G_T =
phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dT;
GlobalDimMatrixType const diffusion_W_L_T =
@@ -1689,25 +1672,19 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
temperature_index)
.noalias() += gradNpT * ip_dd.dfW_4_LWpG.dT * gradpGR * NT * w;
- LWpC.noalias() -=
- gradNpT * (ip_cv.advection_data.advection_W_L + diffusion_W_L_p) *
- gradNp * w;
+ LWpC.noalias() += gradNpT * ip_cv.fW_4_LWpC.L * gradNp * w;
// fW_4 LWp_cap' parts; LWpC = \int grad (a + d) grad
local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
- gradNpT * (ip_cv.dfW_4_LWpC_a_dp_GR + ip_cv.dfW_4_LWpC_d_dp_GR) *
- gradpCap * Np * w;
+ gradNpT * ip_dd.dfW_4_LWpC.dp_GR * gradpCap * Np * w;
local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
- gradNpT * (ip_cv.dfW_4_LWpC_a_dp_cap + ip_cv.dfW_4_LWpC_d_dp_cap) *
- gradpCap * Np * w;
+ gradNpT * ip_dd.dfW_4_LWpC.dp_cap * gradpCap * Np * w;
local_Jac
.template block<W_size, temperature_size>(W_index,
temperature_index)
- .noalias() -= gradNpT *
- (ip_cv.dfW_4_LWpC_a_dT + ip_cv.dfW_4_LWpC_d_dT) *
- gradpCap * NT * w;
+ .noalias() -= gradNpT * ip_dd.dfW_4_LWpC.dT * gradpCap * NT * w;
LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;