From a56e19a541820dd4973977bf0a816afc2b79d14b Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Wed, 3 Apr 2024 18:50:29 +0200 Subject: [PATCH 01/22] [PL/TH2M] Extract solid heat capacity data/model --- .../ConstitutiveRelations/ConstitutiveData.h | 2 ++ .../ConstitutiveModels.h | 2 ++ .../SolidHeatCapacity.cpp | 36 +++++++++++++++++++ .../ConstitutiveRelations/SolidHeatCapacity.h | 32 +++++++++++++++++ ProcessLib/TH2M/TH2MFEM-impl.h | 14 ++++---- 5 files changed, 80 insertions(+), 6 deletions(-) create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/SolidHeatCapacity.cpp create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/SolidHeatCapacity.h diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h index 0965f8009c6..f79ba86ea6b 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h @@ -28,6 +28,7 @@ #include "Saturation.h" #include "SolidCompressibility.h" #include "SolidDensity.h" +#include "SolidHeatCapacity.h" #include "SolidMechanics.h" #include "SolidThermalExpansion.h" #include "Swelling.h" @@ -138,6 +139,7 @@ struct ConstitutiveTempData PhaseTransitionData phase_transition_data; PorosityDerivativeData porosity_d_data; SolidDensityDerivativeData solid_density_d_data; + SolidHeatCapacityData solid_heat_capacity_data; using DisplacementDimVector = Eigen::Matrix; using DisplacementDimMatrix = diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h index b9526cdf6d3..50fedd23598 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h @@ -20,6 +20,7 @@ #include "Saturation.h" #include "SolidCompressibility.h" #include "SolidDensity.h" +#include "SolidHeatCapacity.h" #include "SolidMechanics.h" #include "SolidThermalExpansion.h" #include "Swelling.h" @@ -69,6 +70,7 @@ struct ConstitutiveModels PorosityModel porosity_model; SolidDensityModel solid_density_model; #endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION + SolidHeatCapacityModel solid_heat_capacity_model; }; } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidHeatCapacity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/SolidHeatCapacity.cpp new file mode 100644 index 00000000000..7e1fcef50cd --- /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(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 00000000000..4f3b4fae0ad --- /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; + +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/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h index 9203a3d5326..df952164fb7 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -246,6 +246,9 @@ TH2MLocalAssembler(vars, pos, t, dt); - ip_data.h_S = cpS * T; + ip_data.h_S = ip_cv.solid_heat_capacity_data() * 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 + @@ -348,7 +348,8 @@ TH2MLocalAssembler Date: Wed, 3 Apr 2024 19:31:10 +0200 Subject: [PATCH 02/22] [PL/TH2M] Collect thermal conductivities data --- ...fectiveThermalConductivityPorosityMixing.h | 6 ++-- .../ConstitutiveRelations/ConstitutiveData.h | 5 ++- .../ThermalConductivity.h | 31 ++++++++++++++++++ ProcessLib/TH2M/IntegrationPointData.h | 1 - ProcessLib/TH2M/TH2MFEM-impl.h | 32 +++++++++++-------- 5 files changed, 55 insertions(+), 20 deletions(-) create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h diff --git a/MaterialLib/MPL/Properties/EffectiveThermalConductivityPorosityMixing.h b/MaterialLib/MPL/Properties/EffectiveThermalConductivityPorosityMixing.h index f5ef7fdbd78..01f80c95962 100644 --- a/MaterialLib/MPL/Properties/EffectiveThermalConductivityPorosityMixing.h +++ b/MaterialLib/MPL/Properties/EffectiveThermalConductivityPorosityMixing.h @@ -39,9 +39,9 @@ class EffectiveThermalConductivityPorosityMixing final : public Property 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/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h index f79ba86ea6b..a543613ec98 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h @@ -32,6 +32,7 @@ #include "SolidMechanics.h" #include "SolidThermalExpansion.h" #include "Swelling.h" +#include "ThermalConductivity.h" #include "TotalStress.h" #include "VapourPartialPressure.h" #include "Viscosity.h" @@ -140,14 +141,12 @@ struct ConstitutiveTempData PorosityDerivativeData porosity_d_data; SolidDensityDerivativeData solid_density_d_data; SolidHeatCapacityData solid_heat_capacity_data; + ThermalConductivityData thermal_conductivity_data; using DisplacementDimVector = Eigen::Matrix; using DisplacementDimMatrix = Eigen::Matrix; - 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; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h new file mode 100644 index 00000000000..2ec711957f0 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h @@ -0,0 +1,31 @@ +/** + * \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" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ + +template +struct ThermalConductivityData +{ + GlobalDimMatrix lambda; + // Currently unused, but there is a comment in TH2MFEM-impl.h referring to + // this matrix + // GlobalDimMatrix dlambda_dp_GR; + GlobalDimMatrix dlambda_dp_cap; + GlobalDimMatrix dlambda_dT; +}; + +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h index ec9de5fdb91..e33c29163ed 100644 --- a/ProcessLib/TH2M/IntegrationPointData.h +++ b/ProcessLib/TH2M/IntegrationPointData.h @@ -48,7 +48,6 @@ struct IntegrationPointData final double rho_u_eff = std::numeric_limits::quiet_NaN(); double rho_u_eff_prev = std::numeric_limits::quiet_NaN(); - GlobalDimMatrixType lambda; GlobalDimVectorType d_CG; GlobalDimVectorType d_WG; GlobalDimVectorType d_CL; diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h index df952164fb7..380f15766ee 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -273,11 +273,12 @@ TH2MLocalAssembler( - medium - .property( - MaterialPropertyLib::PropertyType::thermal_conductivity) - .value(vars, pos, t, dt)); + ip_cv.thermal_conductivity_data.lambda = + MaterialPropertyLib::formEigenTensor( + medium + .property( + MaterialPropertyLib::PropertyType::thermal_conductivity) + .value(vars, pos, t, dt)); ip_data.h_S = ip_cv.solid_heat_capacity_data() * T; auto const u_S = ip_data.h_S; @@ -409,12 +410,12 @@ TH2MLocalAssembler(0.); - ip_cv.dlambda_dp_cap = + ip_cv.thermal_conductivity_data.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; + ip_cv.thermal_conductivity_data.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 @@ -1392,7 +1393,8 @@ void TH2MLocalAssembler< MTu.noalias() += NTT * rho_h_eff * mT * Bu * w; - KTT.noalias() += gradNTT * ip.lambda * gradNT * w; + KTT.noalias() += + gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w; fT.noalias() -= NTT * rho_u_eff_dot * w; @@ -2053,7 +2055,8 @@ void TH2MLocalAssembler(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_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 fT.noalias() -= NTT * rho_u_eff_dot * w; From 2b00d20e90e54e6450d5b9d2f8e3c018056e3d62 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Wed, 3 Apr 2024 20:13:32 +0200 Subject: [PATCH 03/22] [PL/TH2M] Extract thermal cond. computations Both, the value and the derivatives. --- .../ConstitutiveModels.h | 2 + .../ThermalConductivity.cpp | 105 ++++++++++++++++++ .../ThermalConductivity.h | 16 +++ ProcessLib/TH2M/TH2MFEM-impl.h | 67 +---------- 4 files changed, 128 insertions(+), 62 deletions(-) create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h index 50fedd23598..ec38e07b1a9 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h @@ -24,6 +24,7 @@ #include "SolidMechanics.h" #include "SolidThermalExpansion.h" #include "Swelling.h" +#include "ThermalConductivity.h" #include "TotalStress.h" #include "Viscosity.h" @@ -71,6 +72,7 @@ struct ConstitutiveModels SolidDensityModel solid_density_model; #endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION SolidHeatCapacityModel solid_heat_capacity_model; + ThermalConductivityModel thermal_conductivity_model; }; } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp new file mode 100644 index 00000000000..ef443a87ff6 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp @@ -0,0 +1,105 @@ +/** + * \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 +void ThermalConductivityModel::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& 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( + 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( + media_data.gas[MPL::PropertyType::thermal_conductivity].value( + variables, x_t.x, x_t.t, x_t.dt)) + : MPL::formEigenTensor(0.); + + auto const dlambda_GR_dT = + media_data.gas.hasProperty(MPL::PropertyType::thermal_conductivity) + ? MPL::formEigenTensor( + media_data.gas[MPL::PropertyType::thermal_conductivity] + .dValue(variables, MPL::Variable::temperature, x_t.x, + x_t.t, x_t.dt)) + : MPL::formEigenTensor(0.); + + auto const lambdaLR = + media_data.liquid.hasProperty(MPL::PropertyType::thermal_conductivity) + ? MPL::formEigenTensor( + media_data.liquid[MPL::PropertyType::thermal_conductivity] + .value(variables, x_t.x, x_t.t, x_t.dt)) + : MPL::formEigenTensor(0.); + + auto const dlambda_LR_dT = + media_data.liquid.hasProperty(MPL::PropertyType::thermal_conductivity) + ? MPL::formEigenTensor( + media_data.liquid[MPL::PropertyType::thermal_conductivity] + .dValue(variables, MPL::Variable::temperature, x_t.x, + x_t.t, x_t.dt)) + : MPL::formEigenTensor(0.); + + auto const lambdaSR = + media_data.solid.hasProperty(MPL::PropertyType::thermal_conductivity) + ? MPL::formEigenTensor( + media_data.solid[MPL::PropertyType::thermal_conductivity] + .value(variables, x_t.x, x_t.t, x_t.dt)) + : MPL::formEigenTensor(0.); + + auto const dlambda_SR_dT = + media_data.solid.hasProperty(MPL::PropertyType::thermal_conductivity) + ? MPL::formEigenTensor( + media_data.solid[MPL::PropertyType::thermal_conductivity] + .dValue(variables, MPL::Variable::temperature, x_t.x, + x_t.t, x_t.dt)) + : MPL::formEigenTensor(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; + + thermal_conductivity_data.dlambda_dT = + phi_G * dlambda_GR_dT + phi_L * dlambda_LR_dT + phi_S * dlambda_SR_dT - + porosity_d_data.dphi_dT * lambdaSR; +} +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 index 2ec711957f0..fbf34021e55 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h @@ -10,6 +10,8 @@ #pragma once #include "Base.h" +#include "Porosity.h" +#include "Saturation.h" namespace ProcessLib::TH2M { @@ -27,5 +29,19 @@ struct ThermalConductivityData GlobalDimMatrix dlambda_dT; }; +template +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& + thermal_conductivity_data) const; +}; + +extern template struct ThermalConductivityModel<2>; +extern template struct ThermalConductivityModel<3>; } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h index 380f15766ee..83835eb5dc8 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -115,9 +115,7 @@ TH2MLocalAssemblerprocess_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 = @@ -249,6 +247,11 @@ TH2MLocalAssembler( - medium - .property( - MaterialPropertyLib::PropertyType::thermal_conductivity) - .value(vars, pos, t, dt)); - ip_data.h_S = ip_cv.solid_heat_capacity_data() * T; auto const u_S = ip_data.h_S; @@ -365,58 +360,6 @@ TH2MLocalAssembler( - gas_phase - .property(MPL::PropertyType::thermal_conductivity) - .value(vars, pos, t, dt)) - : MPL::formEigenTensor(0.); - - auto const dlambda_GR_dT = - gas_phase.hasProperty(MPL::PropertyType::thermal_conductivity) - ? MPL::formEigenTensor( - gas_phase[MPL::PropertyType::thermal_conductivity].dValue( - vars, MPL::Variable::temperature, pos, t, dt)) - : MPL::formEigenTensor(0.); - - auto const lambdaLR = - liquid_phase.hasProperty(MPL::PropertyType::thermal_conductivity) - ? MPL::formEigenTensor( - liquid_phase - .property(MPL::PropertyType::thermal_conductivity) - .value(vars, pos, t, dt)) - : MPL::formEigenTensor(0.); - - auto const dlambda_LR_dT = - liquid_phase.hasProperty(MPL::PropertyType::thermal_conductivity) - ? MPL::formEigenTensor( - liquid_phase[MPL::PropertyType::thermal_conductivity] - .dValue(vars, MPL::Variable::temperature, pos, t, dt)) - : MPL::formEigenTensor(0.); - - auto const lambdaSR = - solid_phase.hasProperty(MPL::PropertyType::thermal_conductivity) - ? MPL::formEigenTensor( - solid_phase - .property(MPL::PropertyType::thermal_conductivity) - .value(vars, pos, t, dt)) - : MPL::formEigenTensor(0.); - - auto const dlambda_SR_dT = - solid_phase.hasProperty(MPL::PropertyType::thermal_conductivity) - ? MPL::formEigenTensor( - solid_phase[MPL::PropertyType::thermal_conductivity] - .dValue(vars, MPL::Variable::temperature, pos, t, dt)) - : MPL::formEigenTensor(0.); - - ip_cv.thermal_conductivity_data.dlambda_dp_cap = - dphi_G_dp_cap * lambdaGR + dphi_L_dp_cap * lambdaLR; - - ip_cv.thermal_conductivity_data.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) From 8170fb10d893039f4aa87a057a54558dc018d3a7 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Thu, 4 Apr 2024 11:34:39 +0200 Subject: [PATCH 04/22] [PL/TH2M] Default initialize scalars to NaN. --- ProcessLib/TH2M/ConstitutiveRelations/Base.h | 4 ++-- ProcessLib/TH2M/ConstitutiveRelations/FluidDensity.h | 4 ++-- ProcessLib/TH2M/ConstitutiveRelations/MassMoleFractions.h | 8 ++++---- .../TH2M/ConstitutiveRelations/SolidThermalExpansion.h | 4 ++-- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Base.h b/ProcessLib/TH2M/ConstitutiveRelations/Base.h index d0f26028eca..a2bcea4a0dc 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/Base.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/Base.h @@ -47,8 +47,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/FluidDensity.h b/ProcessLib/TH2M/ConstitutiveRelations/FluidDensity.h index 067bc887273..22a0fa60d87 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/MassMoleFractions.h b/ProcessLib/TH2M/ConstitutiveRelations/MassMoleFractions.h index a63b2514a47..a180d179d83 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/SolidThermalExpansion.h b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h index 9f225f30f7d..3eeab68c3d2 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h @@ -19,8 +19,8 @@ template struct SolidThermalExpansionData { KelvinVector solid_linear_thermal_expansivity_vector; - double beta_T_SR; /// Isotropic solid phase volumetric thermal expansion - /// coefficient. + double beta_T_SR = nan; /// Isotropic solid phase volumetric thermal + /// expansion coefficient. double thermal_volume_strain = nan; }; From 960026b1fee42f53a51e254879fa36191efa688c Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Thu, 25 Apr 2024 18:48:50 +0200 Subject: [PATCH 05/22] [PL/TH2M] Improve comment --- ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h index 3eeab68c3d2..6d4a69a73d1 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h @@ -19,8 +19,8 @@ template struct SolidThermalExpansionData { KelvinVector solid_linear_thermal_expansivity_vector; - double beta_T_SR = nan; /// Isotropic solid phase volumetric thermal - /// expansion coefficient. + double beta_T_SR = nan; /// Solid phase volumetric thermal expansion + /// coefficient. double thermal_volume_strain = nan; }; From b6f24b1cd88588c5ee0cac63ab1daf8d2836f7bf Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Thu, 4 Apr 2024 13:05:01 +0200 Subject: [PATCH 06/22] [PL/TH2M] Use drho/dT derivative from phase trans. The computation of the derivative is different in both phase transition models. This removes the last part of MPL evaluation. --- .../NoPhaseTransition.cpp | 4 +-- .../ConstitutiveRelations/PermeabilityData.h | 8 +++--- .../ConstitutiveRelations/PhaseTransition.cpp | 4 +-- ProcessLib/TH2M/TH2MFEM-impl.h | 28 ++----------------- 4 files changed, 11 insertions(+), 33 deletions(-) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp index 9c01e89ce41..974c4c24ec3 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp @@ -152,7 +152,7 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, cv.drho_C_LR_dp_GR = 0; cv.drho_C_GR_dT = drho_GR_dT; - auto const drho_LR_dT = + cv.drho_LR_dT = liquid_phase[MaterialPropertyLib::PropertyType::density] .template dValue(variables, MaterialPropertyLib::Variable::temperature, @@ -169,7 +169,7 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, 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_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; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PermeabilityData.h b/ProcessLib/TH2M/ConstitutiveRelations/PermeabilityData.h index b6eeac5c148..aea68ad3d26 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PermeabilityData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/PermeabilityData.h @@ -19,10 +19,10 @@ namespace ConstitutiveRelations template 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 Ki; static auto reflect() diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp index b34faedf641..48d922cfaaf 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp @@ -465,7 +465,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 +480,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 = diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h index 83835eb5dc8..e7e6336c4d7 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -115,7 +115,6 @@ TH2MLocalAssemblerprocess_data_.media_map.getMedium(this->element_.getID()); - auto const& liquid_phase = medium.phase("AqueousLiquid"); ConstitutiveRelations::MediaData media_data{medium}; unsigned const n_integration_points = @@ -160,7 +159,6 @@ TH2MLocalAssemblerprocess_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; @@ -252,21 +250,6 @@ TH2MLocalAssemblerS_L; - - vars.porosity = ip_out.porosity_data.phi; - auto const& c = ip_cv.phase_transition_data; auto const phi_L = @@ -333,15 +316,10 @@ TH2MLocalAssembler(vars, MPL::Variable::temperature, pos, - t, dt); - ip_cv.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 * @@ -393,7 +371,7 @@ TH2MLocalAssembler Date: Thu, 4 Apr 2024 16:14:32 +0200 Subject: [PATCH 07/22] [PL/TH2M] Make hCL and hWL local to PhaseTrans. Since both variables are used only internally one can same a little space here. --- ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp | 6 +++--- ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h | 2 -- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp index 48d922cfaaf..f809791cd05 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp @@ -508,9 +508,9 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, .template value(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; // specific inner energies of liquid phase cv.uL = enthalpy_data.h_L; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h index 2500aa882ae..f1e14c949e2 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h @@ -64,8 +64,6 @@ struct PhaseTransitionData // specific enthalpies double hCG = 0; double hWG = 0; - double hCL = 0; - double hWL = 0; double dh_G_dT = 0; double dh_L_dT = 0; From 2b22ad7920219b2556f8cd803826b1d006eaa0fd Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Thu, 4 Apr 2024 16:16:29 +0200 Subject: [PATCH 08/22] [PL/TH2M] Store only xnCG avoiding xnWG storage The xnWG is not used in assembly and is required only internally (and in the unit test) but is easily substituted with 1-xnCG. --- .../TH2M/ConstitutiveRelations/NoPhaseTransition.cpp | 3 +-- .../TH2M/ConstitutiveRelations/PhaseTransition.cpp | 8 ++++---- .../TH2M/ConstitutiveRelations/PhaseTransitionData.h | 1 - Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp | 10 ++++++---- 4 files changed, 11 insertions(+), 11 deletions(-) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp index 974c4c24ec3..b9c18f91284 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp @@ -66,9 +66,8 @@ 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.xnCG = 1.; mass_mole_fractions_data.xmCG = 1. - cv.xmWG; auto const M = diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp index f809791cd05..1ea9351d52b 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp @@ -233,9 +233,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,11 +243,11 @@ 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; + cv.xmWG = xnWG * M_W / MG; mass_mole_fractions_data.xmCG = 1. - cv.xmWG; auto const dxn_dxm_conversion = M_W * M_C / MG / MG; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h index f1e14c949e2..5cd94108576 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h @@ -45,7 +45,6 @@ struct PhaseTransitionData double drho_W_LR_dp_LR = 0.; // constituent mass and molar fractions - double xnWG = 0.; double xmWG = 0.; // mass fraction derivatives diff --git a/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp b/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp index 7fbe098c97d..e4a4ed4d45e 100644 --- a/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp +++ b/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp @@ -397,7 +397,8 @@ 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 @@ -424,8 +425,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); @@ -655,7 +656,8 @@ 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 From f26169d49455e70301179740255bd74e2e6edaaf Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Thu, 4 Apr 2024 16:27:14 +0200 Subject: [PATCH 09/22] [PL/TH2M] Storing drho_LR_dp_GR is not neccessary as it is used only locally in the NoPhaseTransition model --- ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp | 3 +-- ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h | 1 - 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp index b9c18f91284..6480d8f6909 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp @@ -140,7 +140,6 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, .template dValue( 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 / @@ -167,7 +166,7 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, */ 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_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; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h index 5cd94108576..8c6a3626189 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h @@ -29,7 +29,6 @@ struct PhaseTransitionData double drho_W_GR_dp_cap = 0.; double drho_W_GR_dT = 0.; - double drho_LR_dp_GR = 0.; double drho_LR_dp_cap = 0.; double drho_LR_dT = 0.; double drho_LR_dp_LR = 0.; From 4bba06a81c8c31f7198a710607c5bc80910a5178 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Thu, 4 Apr 2024 16:28:03 +0200 Subject: [PATCH 10/22] [PL/TH2M] Remove two unused derivatives from PTD PTD = PhaseTransitionData. --- ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h index 8c6a3626189..4680f66bcd8 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h @@ -34,12 +34,10 @@ struct PhaseTransitionData double drho_LR_dp_LR = 0.; 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.; 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.; From 622bcd0bc98fa5720d9f7201273a1446d5c1a184 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Tue, 9 Apr 2024 18:42:48 +0200 Subject: [PATCH 11/22] [PL/TH2M] Remove duplicate unused derivatives Already stored in phase transition data and computed by the corresponding models. --- ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h index a543613ec98..00e55b8cbab 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h @@ -181,8 +181,6 @@ struct ConstitutiveTempData double drho_h_eff_dT = std::numeric_limits::quiet_NaN(); double drho_h_eff_dp_GR = std::numeric_limits::quiet_NaN(); double drho_h_eff_dp_cap = std::numeric_limits::quiet_NaN(); - double dh_G_dT = std::numeric_limits::quiet_NaN(); - double dh_L_dT = std::numeric_limits::quiet_NaN(); double dfC_4_MCpG_dp_GR = std::numeric_limits::quiet_NaN(); double dfC_4_MCpG_dT = std::numeric_limits::quiet_NaN(); double dfC_4_MCT_dT = std::numeric_limits::quiet_NaN(); From e67e3db450531d9f869fadf0318bca5c87760ef3 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Thu, 4 Apr 2024 17:51:17 +0200 Subject: [PATCH 12/22] [PL/TH2M] NaN initialize PhaseTransitionData - Some variables are set in one model but not in the other: these default 0 initializations are now explicit in both models. - Some variables are not initialized or initialized to zero in both models: these are marked as constexpr 0 in the PhaseTransitionData with a TODO note. --- .../NoPhaseTransition.cpp | 32 ++++--- .../ConstitutiveRelations/PhaseTransition.cpp | 11 +++ .../PhaseTransitionData.h | 89 ++++++++++--------- 3 files changed, 78 insertions(+), 54 deletions(-) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp index 6480d8f6909..3a98aeab783 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp @@ -67,6 +67,9 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, // C-component is only component in the gas phase cv.xmWG = 0.; + cv.dxmWG_dpGR = 0.; + cv.dxmWG_dpCap = 0.; + cv.dxmWG_dT = 0.; mass_mole_fractions_data.xnCG = 1.; mass_mole_fractions_data.xmCG = 1. - cv.xmWG; @@ -89,6 +92,9 @@ 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; @@ -120,16 +126,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(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( @@ -146,24 +154,14 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, 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; cv.drho_LR_dT = liquid_phase[MaterialPropertyLib::PropertyType::density] .template dValue(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_LR; @@ -171,6 +169,12 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, 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/PhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp index 1ea9351d52b..59901703477 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp @@ -337,9 +337,12 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, // specific enthalpy of gas phase enthalpy_data.h_G = mass_mole_fractions_data.xmCG * cv.hCG + cv.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 = @@ -511,9 +514,11 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, 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 = @@ -526,6 +531,12 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, viscosity_data.mu_LR = liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity) .template value(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/PhaseTransitionData.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h index 4680f66bcd8..98cd682b4ee 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h @@ -17,62 +17,71 @@ 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_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_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_dT = 0.; - double drho_W_LR_dp_LR = 0.; + double drho_W_LR_dp_GR = nan; + double drho_W_LR_dT = nan; + double drho_W_LR_dp_LR = nan; // constituent mass and molar fractions - double xmWG = 0.; + double xmWG = 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 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 From 8526548679ef1c93c724c5fb6f04d2601632d0a8 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Thu, 4 Apr 2024 18:36:14 +0200 Subject: [PATCH 13/22] [PL/TH2M] Replace xmWG with 1 - xmCG Reducing storage a little more. --- .../NoPhaseTransition.cpp | 3 +- .../ConstitutiveRelations/PhaseTransition.cpp | 17 +++++----- .../PhaseTransitionData.h | 3 -- ProcessLib/TH2M/TH2MFEM-impl.h | 11 ++++--- .../TH2M/TestTH2MNoPhaseTransition.cpp | 3 +- .../TH2M/TestTH2MPhaseTransition.cpp | 32 +++++++++---------- 6 files changed, 32 insertions(+), 37 deletions(-) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp index 3a98aeab783..43f2212b006 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp @@ -66,12 +66,11 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, vapour_pressure_data.pWGR = 0; // C-component is only component in the gas phase - cv.xmWG = 0.; cv.dxmWG_dpGR = 0.; cv.dxmWG_dpCap = 0.; cv.dxmWG_dT = 0.; mass_mole_fractions_data.xnCG = 1.; - mass_mole_fractions_data.xmCG = 1. - cv.xmWG; + mass_mole_fractions_data.xmCG = 1.; auto const M = gas_phase.property(MaterialPropertyLib::PropertyType::molar_mass) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp index 59901703477..748fd9cca98 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp @@ -247,8 +247,8 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, variables.molar_mass = MG; // gas phase mass fractions - cv.xmWG = 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 +302,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 +314,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,8 +335,7 @@ 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 diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h index 98cd682b4ee..a0042c6afbb 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h @@ -43,9 +43,6 @@ struct PhaseTransitionData double drho_W_LR_dT = nan; double drho_W_LR_dp_LR = nan; - // constituent mass and molar fractions - double xmWG = nan; - // mass fraction derivatives double dxmWG_dpGR = nan; double dxmWG_dpCap = nan; diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h index e7e6336c4d7..edad8484191 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -295,11 +295,12 @@ TH2MLocalAssembler Date: Fri, 5 Apr 2024 14:12:46 +0200 Subject: [PATCH 14/22] [PL/TH2M] Extract viscosity model from phase model From both phase transition models. --- .../ConstitutiveModels.h | 1 + .../NoPhaseTransition.cpp | 8 -- .../ConstitutiveRelations/NoPhaseTransition.h | 2 +- .../ConstitutiveRelations/PhaseTransition.cpp | 16 +--- .../ConstitutiveRelations/PhaseTransition.h | 2 +- .../PhaseTransitionModel.h | 6 +- .../TH2M/ConstitutiveRelations/Viscosity.cpp | 38 ++++++++ .../TH2M/ConstitutiveRelations/Viscosity.h | 9 ++ ProcessLib/TH2M/TH2MFEM-impl.h | 6 +- .../TH2M/TestTH2MNoPhaseTransition.cpp | 14 +-- .../TH2M/TestTH2MPhaseTransition.cpp | 88 ++++++++----------- 11 files changed, 97 insertions(+), 93 deletions(-) create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/Viscosity.cpp diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h index ec38e07b1a9..e37d7c906d2 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h @@ -62,6 +62,7 @@ struct ConstitutiveModels PermeabilityModel permeability_model; PureLiquidDensityModel pure_liquid_density_model; PhaseTransitionModel const& phase_transition_model; + ViscosityModel viscosity_model; #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION PorosityModelNonConstantSolidPhaseVolumeFraction porosity_model; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp index 43f2212b006..782b8b60ecb 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, @@ -81,9 +80,6 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, fluid_density_data.rho_GR = gas_phase.property(MaterialPropertyLib::PropertyType::density) .template value(variables, x_t.x, x_t.t, x_t.dt); - viscosity_data.mu_GR = - gas_phase.property(MaterialPropertyLib::PropertyType::viscosity) - .template value(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; @@ -100,10 +96,6 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, 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(variables, x_t.x, x_t.t, x_t.dt); - // specific heat capacities auto const cpG = gas_phase diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h index 20c7b5abd0a..918723e9060 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/PhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp index 748fd9cca98..7f71fcb2492 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, @@ -357,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(variables, x_t.x, x_t.t, x_t.dt); - // Dissolution part -- Liquid phase properties // ------------------------------------------- @@ -526,11 +519,6 @@ 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(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.; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h index d4ffbc950f2..69d621b3dad 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/PhaseTransitionModel.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h index 63672d93664..8cf92e94f7e 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/Viscosity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Viscosity.cpp new file mode 100644 index 00000000000..782dcb40c5c --- /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(variables, x_t.x, x_t.t, x_t.dt); + + viscosity_data.mu_LR = + liquid_phase[MaterialPropertyLib::PropertyType::viscosity] + .template value(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 95a38fd1d24..c46e48c702d 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/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h index edad8484191..1ecdc9cfd6d 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -223,11 +223,15 @@ TH2MLocalAssembler\n"; m << " \n"; @@ -54,7 +51,6 @@ TEST(ProcessLib, TH2MNoPhaseTransition) // gas phase properties m << "\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 << "\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,8 +135,6 @@ 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); @@ -163,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 3bdf25f935f..f636edf49f6 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 << "\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 << "\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,9 +241,8 @@ 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 = 1 - mass_mole_fractions.xmCG; auto rhoGR_plus = fluid_density.rho_GR; @@ -258,9 +252,8 @@ 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 = 1 - mass_mole_fractions.xmCG; auto rhoGR_minus = fluid_density.rho_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,8 +290,8 @@ 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 = 1 - mass_mole_fractions.xmCG; rhoGR_plus = fluid_density.rho_GR; @@ -310,8 +303,8 @@ 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 = 1 - mass_mole_fractions.xmCG; rhoGR_minus = fluid_density.rho_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,8 +341,8 @@ 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 = 1 - mass_mole_fractions.xmCG; @@ -360,8 +353,8 @@ 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 = 1 - mass_mole_fractions.xmCG; @@ -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); @@ -450,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); } } @@ -490,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; @@ -514,9 +503,8 @@ 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 = 1 - mass_mole_fractions.xmCG; auto rhoCGR_plus = constituent_density.rho_C_GR; @@ -525,9 +513,8 @@ 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_minus = 1 - mass_mole_fractions.xmCG; auto rhoCGR_minus = constituent_density.rho_C_GR; @@ -537,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); @@ -562,8 +549,8 @@ 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 = 1 - mass_mole_fractions.xmCG; rhoCGR_plus = constituent_density.rho_C_GR; @@ -574,8 +561,8 @@ 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 = 1 - mass_mole_fractions.xmCG; rhoCGR_minus = constituent_density.rho_C_GR; @@ -585,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); @@ -610,8 +597,8 @@ 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 = 1 - mass_mole_fractions.xmCG; @@ -621,8 +608,8 @@ 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_minus = 1 - mass_mole_fractions.xmCG; @@ -633,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); @@ -706,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); } } From 167f3de86b159646a3478298ecde214483f76628 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 5 Apr 2024 12:52:20 +0200 Subject: [PATCH 15/22] [PL/TH2M] Move h_S together with other enthalpies --- .../TH2M/ConstitutiveRelations/Enthalpy.h | 4 +++- ProcessLib/TH2M/IntegrationPointData.h | 3 --- ProcessLib/TH2M/LocalAssemblerInterface.h | 6 ------ ProcessLib/TH2M/TH2MFEM-impl.h | 11 ++++++----- ProcessLib/TH2M/TH2MFEM.h | 18 ++++-------------- ProcessLib/TH2M/TH2MProcess.cpp | 4 ---- 6 files changed, 13 insertions(+), 33 deletions(-) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h index 48dbbaeb764..fa8a4db5fb6 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,7 +28,8 @@ 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)}; } }; diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h index e33c29163ed..8e1e7da6631 100644 --- a/ProcessLib/TH2M/IntegrationPointData.h +++ b/ProcessLib/TH2M/IntegrationPointData.h @@ -41,9 +41,6 @@ struct IntegrationPointData final double rho_L_h_L = std::numeric_limits::quiet_NaN(); double rho_S_h_S = std::numeric_limits::quiet_NaN(); - // specific enthalpies - double h_S = std::numeric_limits::quiet_NaN(); - // internal energies double rho_u_eff = std::numeric_limits::quiet_NaN(); double rho_u_eff_prev = std::numeric_limits::quiet_NaN(); diff --git a/ProcessLib/TH2M/LocalAssemblerInterface.h b/ProcessLib/TH2M/LocalAssemblerInterface.h index 9e8d86a0863..cc0331f44b6 100644 --- a/ProcessLib/TH2M/LocalAssemblerInterface.h +++ b/ProcessLib/TH2M/LocalAssemblerInterface.h @@ -102,12 +102,6 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, std::vector const& dof_table, std::vector& cache) const = 0; - virtual std::vector const& getIntPtEnthalpySolid( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& 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 1ecdc9cfd6d..d00cee689e1 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -262,8 +262,8 @@ TH2MLocalAssembler return Eigen::Map(N_u.data(), N_u.size()); } + // TODO: Here is some refactoring potential. All secondary variables could + // be stored in some container to avoid defining one method for each + // variable. + std::vector const& getIntPtDarcyVelocityGas( const double t, std::vector const& x, @@ -231,20 +235,6 @@ class TH2MLocalAssembler : public LocalAssemblerInterface ConstitutiveRelations::ConstitutiveModels 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 const& getIntPtEnthalpySolid( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const override - { - return ProcessLib::getIntegrationPointScalarData(_ip_data, &IpData::h_S, - cache); - } - private: using BMatricesType = BMatrixPolicyType; diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp index 34d2e0d1db4..6eacdc8c27f 100644 --- a/ProcessLib/TH2M/TH2MProcess.cpp +++ b/ProcessLib/TH2M/TH2MProcess.cpp @@ -186,10 +186,6 @@ void TH2MProcess::initializeConcreteProcess( &LocalAssemblerInterface< DisplacementDim>::getIntPtDiffusionVelocityLiquidLiquid); - add_secondary_variable( - "enthalpy_solid", 1, - &LocalAssemblerInterface::getIntPtEnthalpySolid); - // // enable output of internal variables defined by material models // From 9cbdb4f504adf96f2a6a8800d64d52f507c67b2c Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 5 Apr 2024 13:26:42 +0200 Subject: [PATCH 16/22] [PL/TH2M] Move internal energies to states This also eliminates the IntegrationPointData::pushBackState(). --- .../ConstitutiveRelations/ConstitutiveData.h | 4 ++++ .../ConstitutiveRelations/InternalEnergy.h | 22 +++++++++++++++++++ ProcessLib/TH2M/IntegrationPointData.h | 6 ----- ProcessLib/TH2M/TH2MFEM-impl.h | 18 ++++++++------- ProcessLib/TH2M/TH2MFEM.h | 2 -- 5 files changed, 36 insertions(+), 16 deletions(-) create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h index 00e55b8cbab..f56b9186b2a 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h @@ -16,6 +16,7 @@ #include "Enthalpy.h" #include "EquivalentPlasticStrainData.h" #include "FluidDensity.h" +#include "InternalEnergy.h" #include "MassMoleFractions.h" #include "MechanicalStrain.h" #include "PermeabilityData.h" @@ -52,6 +53,7 @@ struct StatefulData MechanicalStrainData mechanical_strain_data; PureLiquidDensityData rho_W_LR; ConstituentDensityData constituent_density_data; + InternalEnergyData internal_energy_data; static auto reflect() { @@ -72,6 +74,7 @@ struct StatefulDataPrev PrevState> mechanical_strain_data; PrevState rho_W_LR; PrevState constituent_density_data; + PrevState internal_energy_data; StatefulDataPrev& operator=( StatefulData const& state) @@ -82,6 +85,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; } diff --git a/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h new file mode 100644 index 00000000000..f5f393704bc --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h @@ -0,0 +1,22 @@ +/** + * \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; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h index 8e1e7da6631..5a5c1b0d9e3 100644 --- a/ProcessLib/TH2M/IntegrationPointData.h +++ b/ProcessLib/TH2M/IntegrationPointData.h @@ -41,10 +41,6 @@ struct IntegrationPointData final double rho_L_h_L = std::numeric_limits::quiet_NaN(); double rho_S_h_S = std::numeric_limits::quiet_NaN(); - // internal energies - double rho_u_eff = std::numeric_limits::quiet_NaN(); - double rho_u_eff_prev = std::numeric_limits::quiet_NaN(); - GlobalDimVectorType d_CG; GlobalDimVectorType d_WG; GlobalDimVectorType d_CL; @@ -55,8 +51,6 @@ struct IntegrationPointData final double integration_weight = std::numeric_limits::quiet_NaN(); - void pushBackState() { rho_u_eff_prev = rho_u_eff; } - EIGEN_MAKE_ALIGNED_OPERATOR_NEW; }; diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h index d00cee689e1..41abf2db569 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -265,9 +265,10 @@ TH2MLocalAssemblermaterial_states_[ip].pushBackState(); this->prev_states_[ip] = this->current_states_[ip]; } @@ -1131,7 +1129,9 @@ void TH2MLocalAssembler< 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; + auto const rho_u_eff_dot = (current_state.internal_energy_data() - + **prev_state.internal_energy_data) / + dt; GlobalDimMatrixType const k_over_mu_G = ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G / @@ -1623,7 +1623,9 @@ void TH2MLocalAssembler this->solid_material_.initializeInternalStateVariables( t, x_position, *material_state.material_state_variables); - ip_data.pushBackState(); material_state.pushBackState(); } @@ -167,7 +166,6 @@ class TH2MLocalAssembler : public LocalAssemblerInterface for (unsigned ip = 0; ip < n_integration_points; ip++) { - _ip_data[ip].pushBackState(); this->material_states_[ip].pushBackState(); } From 42640a0e19cda4241f1562c66f44d32db4af6833 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 5 Apr 2024 14:52:31 +0200 Subject: [PATCH 17/22] [PL/TH2M] Collect rho_h_eff and rho_u_eff and d/dX Data and their derivatives collected in structs. Avoid storage of intermediate phase enthalpies. --- .../ConstitutiveRelations/ConstitutiveData.h | 10 +-- .../TH2M/ConstitutiveRelations/Enthalpy.h | 12 +++ .../ConstitutiveRelations/InternalEnergy.h | 7 ++ ProcessLib/TH2M/IntegrationPointData.h | 5 -- ProcessLib/TH2M/TH2MFEM-impl.h | 75 +++++++++++-------- 5 files changed, 67 insertions(+), 42 deletions(-) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h index f56b9186b2a..6202ea0934c 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h @@ -146,6 +146,10 @@ struct ConstitutiveTempData SolidDensityDerivativeData solid_density_d_data; SolidHeatCapacityData solid_heat_capacity_data; ThermalConductivityData 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; using DisplacementDimMatrix = @@ -179,12 +183,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::quiet_NaN(); - double drho_u_eff_dp_GR = std::numeric_limits::quiet_NaN(); - double drho_u_eff_dp_cap = std::numeric_limits::quiet_NaN(); - double drho_h_eff_dT = std::numeric_limits::quiet_NaN(); - double drho_h_eff_dp_GR = std::numeric_limits::quiet_NaN(); - double drho_h_eff_dp_cap = std::numeric_limits::quiet_NaN(); double dfC_4_MCpG_dp_GR = std::numeric_limits::quiet_NaN(); double dfC_4_MCpG_dT = std::numeric_limits::quiet_NaN(); double dfC_4_MCT_dT = std::numeric_limits::quiet_NaN(); diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h index fa8a4db5fb6..b7a1686618b 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h @@ -33,5 +33,17 @@ struct EnthalpyData } }; +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/InternalEnergy.h b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h index f5f393704bc..73644146b34 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h @@ -18,5 +18,12 @@ namespace ConstitutiveRelations { using InternalEnergyData = BaseLib::StrongType; + +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/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h index 5a5c1b0d9e3..f5376bc7ec2 100644 --- a/ProcessLib/TH2M/IntegrationPointData.h +++ b/ProcessLib/TH2M/IntegrationPointData.h @@ -36,11 +36,6 @@ struct IntegrationPointData final typename ShapeMatricesTypePressure::NodalRowVectorType N_p; typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p; - // phase enthalpies - double rho_G_h_G = std::numeric_limits::quiet_NaN(); - double rho_L_h_L = std::numeric_limits::quiet_NaN(); - double rho_S_h_S = std::numeric_limits::quiet_NaN(); - GlobalDimVectorType d_CG; GlobalDimVectorType d_WG; GlobalDimVectorType d_CL; diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h index 41abf2db569..46cf4514c5d 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -270,11 +270,11 @@ TH2MLocalAssembler(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_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_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_cv.thermal_conductivity_data.lambda * gradNT * w; @@ -2017,26 +2018,38 @@ void TH2MLocalAssembler(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_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_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 * From 2e3e8aa1a40770e8553c05d101ff364e7ba0509f Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 5 Apr 2024 16:08:34 +0200 Subject: [PATCH 18/22] [PL/TH2M] Move diffusion velocities to output data --- ProcessLib/TH2M/ConstitutiveRelations/Base.h | 3 + .../ConstitutiveRelations/ConstitutiveData.h | 5 +- .../ConstitutiveRelations/DiffusionVelocity.h | 42 +++++ ProcessLib/TH2M/IntegrationPointData.h | 5 - ProcessLib/TH2M/LocalAssemblerInterface.h | 24 --- ProcessLib/TH2M/TH2MFEM-impl.h | 156 ++++-------------- ProcessLib/TH2M/TH2MFEM.h | 20 --- ProcessLib/TH2M/TH2MProcess.cpp | 16 -- 8 files changed, 78 insertions(+), 193 deletions(-) create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Base.h b/ProcessLib/TH2M/ConstitutiveRelations/Base.h index a2bcea4a0dc..0667e3730d9 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/Base.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/Base.h @@ -32,6 +32,9 @@ template using GlobalDimMatrix = Eigen::Matrix; +template +using GlobalDimVector = Eigen::Vector; + struct MediaData { explicit MediaData(MaterialPropertyLib::Medium const& medium) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h index 6202ea0934c..847a62eb2d1 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h @@ -12,6 +12,7 @@ #include "Biot.h" #include "Bishops.h" #include "ConstitutiveDensity.h" +#include "DiffusionVelocity.h" #include "ElasticTangentStiffnessData.h" #include "Enthalpy.h" #include "EquivalentPlasticStrainData.h" @@ -103,6 +104,7 @@ struct OutputData VapourPartialPressureData vapour_pressure_data; PorosityData porosity_data; SolidDensityData solid_density_data; + DiffusionVelocityData diffusion_velocity_data; static auto reflect() { @@ -115,7 +117,8 @@ 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); } }; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h new file mode 100644 index 00000000000..25a203ca2c0 --- /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 +struct DiffusionVelocityData +{ + GlobalDimVector d_CG; + GlobalDimVector d_WG; + GlobalDimVector d_CL; + GlobalDimVector d_WL; + + static auto reflect() + { + using Self = DiffusionVelocityData; + 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/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h index f5376bc7ec2..76ff44a96f0 100644 --- a/ProcessLib/TH2M/IntegrationPointData.h +++ b/ProcessLib/TH2M/IntegrationPointData.h @@ -36,11 +36,6 @@ struct IntegrationPointData final typename ShapeMatricesTypePressure::NodalRowVectorType N_p; typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p; - 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 cc0331f44b6..aafc2bddafa 100644 --- a/ProcessLib/TH2M/LocalAssemblerInterface.h +++ b/ProcessLib/TH2M/LocalAssemblerInterface.h @@ -78,30 +78,6 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, std::vector const& dof_table, std::vector& cache) const = 0; - virtual std::vector const& getIntPtDiffusionVelocityVapourGas( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - - virtual std::vector const& getIntPtDiffusionVelocityGasGas( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - - virtual std::vector const& getIntPtDiffusionVelocitySoluteLiquid( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - - virtual std::vector const& getIntPtDiffusionVelocityLiquidLiquid( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& 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 46cf4514c5d..ff308665c59 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -294,30 +294,32 @@ TH2MLocalAssembler const& TH2MLocalAssembler< return cache; } -template -std::vector const& TH2MLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtDiffusionVelocityVapourGas( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - unsigned const n_integration_points = - this->integration_method_.getNumberOfPoints(); - - cache.clear(); - auto cache_matrix = MathLib::createZeroedMatrix>( - 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 -std::vector const& TH2MLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtDiffusionVelocityGasGas( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - unsigned const n_integration_points = - this->integration_method_.getNumberOfPoints(); - - cache.clear(); - auto cache_matrix = MathLib::createZeroedMatrix>( - 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 -std::vector const& TH2MLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtDiffusionVelocitySoluteLiquid( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - unsigned const n_integration_points = - this->integration_method_.getNumberOfPoints(); - - cache.clear(); - auto cache_matrix = MathLib::createZeroedMatrix>( - 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 -std::vector const& TH2MLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtDiffusionVelocityLiquidLiquid( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - unsigned const n_integration_points = - this->integration_method_.getNumberOfPoints(); - - cache.clear(); - auto cache_matrix = MathLib::createZeroedMatrix>( - 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 void TH2MLocalAssembler std::vector const& x, std::vector const& dof_table, std::vector& cache) const override; - std::vector const& getIntPtDiffusionVelocityVapourGas( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - std::vector const& getIntPtDiffusionVelocityGasGas( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - std::vector const& getIntPtDiffusionVelocitySoluteLiquid( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - std::vector const& getIntPtDiffusionVelocityLiquidLiquid( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; std::tuple< std::vector>, diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp index 6eacdc8c27f..51cd4a76d80 100644 --- a/ProcessLib/TH2M/TH2MProcess.cpp +++ b/ProcessLib/TH2M/TH2MProcess.cpp @@ -169,22 +169,6 @@ void TH2MProcess::initializeConcreteProcess( add_secondary_variable( "velocity_liquid", mesh.getDimension(), &LocalAssemblerInterface::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); // // enable output of internal variables defined by material models From d7dd37a6dab427a8e43d5ee36df7c3b95e83e02b Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 5 Apr 2024 16:31:14 +0200 Subject: [PATCH 19/22] [PL/TH2M] Move darcy velocity to output data --- .../ConstitutiveRelations/ConstitutiveData.h | 5 +- .../ConstitutiveRelations/DarcyVelocity.h | 36 ++++++ ProcessLib/TH2M/IntegrationPointData.h | 3 - ProcessLib/TH2M/LocalAssemblerInterface.h | 12 -- ProcessLib/TH2M/TH2MFEM-impl.h | 108 ++++++------------ ProcessLib/TH2M/TH2MFEM.h | 15 --- ProcessLib/TH2M/TH2MProcess.cpp | 7 -- 7 files changed, 75 insertions(+), 111 deletions(-) create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h index 847a62eb2d1..7c4611cea70 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h @@ -12,6 +12,7 @@ #include "Biot.h" #include "Bishops.h" #include "ConstitutiveDensity.h" +#include "DarcyVelocity.h" #include "DiffusionVelocity.h" #include "ElasticTangentStiffnessData.h" #include "Enthalpy.h" @@ -105,6 +106,7 @@ struct OutputData PorosityData porosity_data; SolidDensityData solid_density_data; DiffusionVelocityData diffusion_velocity_data; + DarcyVelocityData darcy_velocity_data; static auto reflect() { @@ -118,7 +120,8 @@ struct OutputData &Self::vapour_pressure_data, &Self::porosity_data, &Self::solid_density_data, - &Self::diffusion_velocity_data); + &Self::diffusion_velocity_data, + &Self::darcy_velocity_data); } }; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h new file mode 100644 index 00000000000..84e66d31d02 --- /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 +struct DarcyVelocityData +{ + GlobalDimVector w_GS; + GlobalDimVector w_LS; + + static auto reflect() + { + using Self = DarcyVelocityData; + 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/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h index 76ff44a96f0..dffc86438d0 100644 --- a/ProcessLib/TH2M/IntegrationPointData.h +++ b/ProcessLib/TH2M/IntegrationPointData.h @@ -36,9 +36,6 @@ struct IntegrationPointData final typename ShapeMatricesTypePressure::NodalRowVectorType N_p; typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p; - GlobalDimVectorType w_GS; - GlobalDimVectorType w_LS; - double integration_weight = std::numeric_limits::quiet_NaN(); EIGEN_MAKE_ALIGNED_OPERATOR_NEW; diff --git a/ProcessLib/TH2M/LocalAssemblerInterface.h b/ProcessLib/TH2M/LocalAssemblerInterface.h index aafc2bddafa..30e624d0b54 100644 --- a/ProcessLib/TH2M/LocalAssemblerInterface.h +++ b/ProcessLib/TH2M/LocalAssemblerInterface.h @@ -66,18 +66,6 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, std::string_view name, double const* values, int const integration_order) = 0; - virtual std::vector const& getIntPtDarcyVelocityGas( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const = 0; - - virtual std::vector const& getIntPtDarcyVelocityLiquid( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& 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 ff308665c59..2f14daaaeeb 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -428,14 +428,17 @@ TH2MLocalAssembler -std::vector const& TH2MLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtDarcyVelocityGas( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - unsigned const n_integration_points = - this->integration_method_.getNumberOfPoints(); - - cache.clear(); - auto cache_matrix = MathLib::createZeroedMatrix>( - 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 -std::vector const& TH2MLocalAssembler< - ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>:: - getIntPtDarcyVelocityLiquid( - const double /*t*/, - std::vector const& /*x*/, - std::vector const& /*dof_table*/, - std::vector& cache) const -{ - unsigned const n_integration_points = - this->integration_method_.getNumberOfPoints(); - - cache.clear(); - auto cache_matrix = MathLib::createZeroedMatrix>( - 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 void TH2MLocalAssembler return Eigen::Map(N_u.data(), N_u.size()); } - // TODO: Here is some refactoring potential. All secondary variables could - // be stored in some container to avoid defining one method for each - // variable. - - std::vector const& getIntPtDarcyVelocityGas( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - std::vector const& getIntPtDarcyVelocityLiquid( - const double t, - std::vector const& x, - std::vector const& dof_table, - std::vector& cache) const override; - std::tuple< std::vector>, std::vector< diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp index 51cd4a76d80..b3289ca9d14 100644 --- a/ProcessLib/TH2M/TH2MProcess.cpp +++ b/ProcessLib/TH2M/TH2MProcess.cpp @@ -163,13 +163,6 @@ void TH2MProcess::initializeConcreteProcess( LocalAssemblerInterface::getReflectionDataForOutput(), _secondary_variables, getExtrapolator(), local_assemblers_); - add_secondary_variable( - "velocity_gas", mesh.getDimension(), - &LocalAssemblerInterface::getIntPtDarcyVelocityGas); - add_secondary_variable( - "velocity_liquid", mesh.getDimension(), - &LocalAssemblerInterface::getIntPtDarcyVelocityLiquid); - // // enable output of internal variables defined by material models // From cfd19ea0ae7620e3a19e7af5537f7026e4835115 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 5 Apr 2024 16:35:27 +0200 Subject: [PATCH 20/22] [PL/TH2M] Cleanup IntegrationPointData Contains now only the shape matrices and the integration weight. ALso remove no longer necessary alignment macros/allocator. --- ProcessLib/TH2M/IntegrationPointData.h | 13 ++----------- ProcessLib/TH2M/TH2MFEM.h | 8 +++----- 2 files changed, 5 insertions(+), 16 deletions(-) diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h index dffc86438d0..0d96e68979e 100644 --- a/ProcessLib/TH2M/IntegrationPointData.h +++ b/ProcessLib/TH2M/IntegrationPointData.h @@ -10,19 +10,12 @@ #pragma once -#include - -#include "MaterialLib/SolidModels/LinearElasticIsotropic.h" -#include "MathLib/KelvinVector.h" -#include "MathLib/LinAlg/Eigen/EigenMapTools.h" -#include "ParameterLib/Parameter.h" - namespace ProcessLib { namespace TH2M { -template +template struct IntegrationPointData final { using GlobalDimMatrixType = @@ -37,8 +30,6 @@ struct IntegrationPointData final typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p; double integration_weight = std::numeric_limits::quiet_NaN(); - - EIGEN_MAKE_ALIGNED_OPERATOR_NEW; }; } // namespace TH2M diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h index ca7e25ccbba..0d88646ae1f 100644 --- a/ProcessLib/TH2M/TH2MFEM.h +++ b/ProcessLib/TH2M/TH2MFEM.h @@ -201,11 +201,9 @@ class TH2MLocalAssembler : public LocalAssemblerInterface private: using BMatricesType = BMatrixPolicyType; - using IpData = - IntegrationPointData; - std::vector> _ip_data; + using IpData = IntegrationPointData; + std::vector _ip_data; SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType> From 235a27b037ac5088a09c5716401acdbb277907a7 Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Fri, 12 Apr 2024 11:55:20 +0200 Subject: [PATCH 21/22] [T/TH2M] Slightly relax ctests' tolerances --- ProcessLib/TH2M/Tests.cmake | 2 +- Tests/Data/TH2M/H2M/OrthotropicSwelling/square.prj | 2 +- .../HM/Confined_Compression/HM_confined_compression_gas.prj | 2 +- .../HM/Confined_Compression/HM_confined_compression_liquid.prj | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/ProcessLib/TH2M/Tests.cmake b/ProcessLib/TH2M/Tests.cmake index 3c737d92e12..cb581c572aa 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 83b487a961f..0726ab51ab7 100644 --- a/Tests/Data/TH2M/H2M/OrthotropicSwelling/square.prj +++ b/Tests/Data/TH2M/H2M/OrthotropicSwelling/square.prj @@ -476,7 +476,7 @@ square_ts_.*.vtu sigma - 2e-15 + 3e-15 0 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 8c95035d44a..91676a26672 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 @@ HM_confined_compression_gas_ts_.*.vtu epsilon - 5e-15 + 6e-15 1e-15 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 698a2dfdd75..97001751570 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 @@ HM_confined_compression_liquid_ts_.*.vtu epsilon - 5e-15 + 6e-15 1e-15 From 39ec3a9ddf0f3c6267845dbbb4f929d3044ee0aa Mon Sep 17 00:00:00 2001 From: Dmitri Naumov Date: Thu, 25 Apr 2024 17:49:51 +0200 Subject: [PATCH 22/22] [PL/TH2M] Correct dlambda/dT derivative --- .../ConstitutiveRelations/ThermalConductivity.cpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp index ef443a87ff6..abd9e0cc9f3 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp @@ -95,9 +95,16 @@ void ThermalConductivityModel::eval( 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 = - phi_G * dlambda_GR_dT + phi_L * dlambda_LR_dT + phi_S * dlambda_SR_dT - - porosity_d_data.dphi_dT * lambdaSR; + (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>;