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/Base.h b/ProcessLib/TH2M/ConstitutiveRelations/Base.h index d0f26028eca..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) @@ -47,8 +50,8 @@ struct MediaData struct TemperatureData { - double T; - double T_prev; + double T = nan; + double T_prev = nan; }; using ReferenceTemperatureData = diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h index 0965f8009c6..7c4611cea70 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h @@ -12,10 +12,13 @@ #include "Biot.h" #include "Bishops.h" #include "ConstitutiveDensity.h" +#include "DarcyVelocity.h" +#include "DiffusionVelocity.h" #include "ElasticTangentStiffnessData.h" #include "Enthalpy.h" #include "EquivalentPlasticStrainData.h" #include "FluidDensity.h" +#include "InternalEnergy.h" #include "MassMoleFractions.h" #include "MechanicalStrain.h" #include "PermeabilityData.h" @@ -28,9 +31,11 @@ #include "Saturation.h" #include "SolidCompressibility.h" #include "SolidDensity.h" +#include "SolidHeatCapacity.h" #include "SolidMechanics.h" #include "SolidThermalExpansion.h" #include "Swelling.h" +#include "ThermalConductivity.h" #include "TotalStress.h" #include "VapourPartialPressure.h" #include "Viscosity.h" @@ -50,6 +55,7 @@ struct StatefulData MechanicalStrainData mechanical_strain_data; PureLiquidDensityData rho_W_LR; ConstituentDensityData constituent_density_data; + InternalEnergyData internal_energy_data; static auto reflect() { @@ -70,6 +76,7 @@ struct StatefulDataPrev PrevState> mechanical_strain_data; PrevState rho_W_LR; PrevState constituent_density_data; + PrevState internal_energy_data; StatefulDataPrev& operator=( StatefulData const& state) @@ -80,6 +87,7 @@ struct StatefulDataPrev mechanical_strain_data = state.mechanical_strain_data; rho_W_LR = state.rho_W_LR; constituent_density_data = state.constituent_density_data; + internal_energy_data = state.internal_energy_data; return *this; } @@ -97,6 +105,8 @@ struct OutputData VapourPartialPressureData vapour_pressure_data; PorosityData porosity_data; SolidDensityData solid_density_data; + DiffusionVelocityData diffusion_velocity_data; + DarcyVelocityData darcy_velocity_data; static auto reflect() { @@ -109,7 +119,9 @@ struct OutputData &Self::fluid_density_data, &Self::vapour_pressure_data, &Self::porosity_data, - &Self::solid_density_data); + &Self::solid_density_data, + &Self::diffusion_velocity_data, + &Self::darcy_velocity_data); } }; @@ -138,14 +150,17 @@ struct ConstitutiveTempData PhaseTransitionData phase_transition_data; PorosityDerivativeData porosity_d_data; SolidDensityDerivativeData solid_density_d_data; + SolidHeatCapacityData solid_heat_capacity_data; + ThermalConductivityData 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 = 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; @@ -174,14 +189,6 @@ struct ConstitutiveTempData DisplacementDimMatrix dadvection_C_dp_cap; DisplacementDimMatrix dk_over_mu_G_dp_cap; DisplacementDimMatrix dk_over_mu_L_dp_cap; - double drho_u_eff_dT = std::numeric_limits::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 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(); diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h index b9526cdf6d3..e37d7c906d2 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h @@ -20,9 +20,11 @@ #include "Saturation.h" #include "SolidCompressibility.h" #include "SolidDensity.h" +#include "SolidHeatCapacity.h" #include "SolidMechanics.h" #include "SolidThermalExpansion.h" #include "Swelling.h" +#include "ThermalConductivity.h" #include "TotalStress.h" #include "Viscosity.h" @@ -60,6 +62,7 @@ struct ConstitutiveModels PermeabilityModel 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; @@ -69,6 +72,8 @@ struct ConstitutiveModels PorosityModel porosity_model; SolidDensityModel solid_density_model; #endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION + SolidHeatCapacityModel solid_heat_capacity_model; + ThermalConductivityModel thermal_conductivity_model; }; } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h new file mode 100644 index 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/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/ConstitutiveRelations/Enthalpy.h b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h index 48dbbaeb764..b7a1686618b 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h @@ -20,6 +20,7 @@ struct EnthalpyData { double h_G = nan; double h_L = nan; + double h_S = nan; static auto reflect() { @@ -27,9 +28,22 @@ struct EnthalpyData namespace R = ProcessLib::Reflection; return std::tuple{R::makeReflectionData("enthalpy_gas", &Self::h_G), - R::makeReflectionData("enthalpy_liquid", &Self::h_L)}; + R::makeReflectionData("enthalpy_liquid", &Self::h_L), + R::makeReflectionData("enthalpy_solid", &Self::h_S)}; } }; +struct EffectiveVolumetricEnthalpy +{ + double rho_h_eff = nan; +}; + +struct EffectiveVolumetricEnthalpyDerivatives +{ + double drho_h_eff_dT = nan; + double drho_h_eff_dp_GR = nan; + double drho_h_eff_dp_cap = nan; +}; + } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/FluidDensity.h b/ProcessLib/TH2M/ConstitutiveRelations/FluidDensity.h index 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/InternalEnergy.h b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h new file mode 100644 index 00000000000..73644146b34 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h @@ -0,0 +1,29 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + */ + +#pragma once + +#include "Base.h" +#include "BaseLib/StrongType.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +using InternalEnergyData = + BaseLib::StrongType; + +struct EffectiveVolumetricInternalEnergyDerivatives +{ + double drho_u_eff_dT = nan; + double drho_u_eff_dp_GR = nan; + double drho_u_eff_dp_cap = nan; +}; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/MassMoleFractions.h b/ProcessLib/TH2M/ConstitutiveRelations/MassMoleFractions.h index 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/NoPhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp index 9c01e89ce41..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, @@ -66,10 +65,11 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, vapour_pressure_data.pWGR = 0; // C-component is only component in the gas phase - cv.xnWG = 0.; - cv.xmWG = 0.; - mass_mole_fractions_data.xnCG = 1. - cv.xnWG; - mass_mole_fractions_data.xmCG = 1. - cv.xmWG; + cv.dxmWG_dpGR = 0.; + cv.dxmWG_dpCap = 0.; + cv.dxmWG_dT = 0.; + mass_mole_fractions_data.xnCG = 1.; + mass_mole_fractions_data.xmCG = 1.; auto const M = gas_phase.property(MaterialPropertyLib::PropertyType::molar_mass) @@ -80,9 +80,6 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, fluid_density_data.rho_GR = gas_phase.property(MaterialPropertyLib::PropertyType::density) .template value(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; @@ -90,16 +87,15 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, // W-component is only component in the liquid phase mass_mole_fractions_data.xmWL = 1.; + cv.dxmWL_dpGR = 0; + cv.dxmWL_dpCap = 0; + cv.dxmWL_dT = 0; auto const pLR = pGR - pCap; variables.liquid_phase_pressure = pLR; fluid_density_data.rho_LR = rho_W_LR(); variables.density = fluid_density_data.rho_LR; - viscosity_data.mu_LR = - liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity) - .template value(variables, x_t.x, x_t.t, x_t.dt); - // specific heat capacities auto const cpG = gas_phase @@ -121,16 +117,18 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, cv.uG = enthalpy_data.h_G - pGR / fluid_density_data.rho_GR; cv.uL = enthalpy_data.h_L; - auto const drho_GR_dT = + cv.drho_GR_dT = gas_phase[MaterialPropertyLib::PropertyType::density] .template dValue(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( @@ -141,38 +139,33 @@ 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 / fluid_density_data.rho_GR; cv.drho_C_GR_dp_GR = cv.drho_GR_dp_GR; - cv.drho_C_LR_dp_LR = 0; - cv.drho_C_LR_dp_GR = 0; - cv.drho_C_GR_dT = drho_GR_dT; + cv.drho_C_GR_dp_cap = 0; + cv.drho_C_GR_dT = cv.drho_GR_dT; - auto const drho_LR_dT = + cv.drho_LR_dT = liquid_phase[MaterialPropertyLib::PropertyType::density] .template dValue(variables, MaterialPropertyLib::Variable::temperature, x_t.x, x_t.t, x_t.dt); - cv.drho_C_LR_dT = 0; - - cv.du_L_dp_GR = 0; - cv.du_L_dp_cap = 0; - /* TODO update to the following when uL has same structure as the uG: - +-1 / fluid_density_data.rho_LR + pLR* cv.drho_LR_dp_cap / - fluid_density_data.rho_LR / - fluid_density_data.rho_LR; - */ cv.drho_W_LR_dp_LR = cv.drho_LR_dp_LR; - cv.drho_W_LR_dp_GR = cv.drho_LR_dp_GR; - cv.drho_W_LR_dT = drho_LR_dT; + cv.drho_W_LR_dp_GR = cv.drho_LR_dp_LR; + cv.drho_W_LR_dT = cv.drho_LR_dT; cv.drho_W_GR_dT = 0; cv.drho_W_GR_dp_GR = 0; cv.drho_W_GR_dp_cap = 0; + + cv.diffusion_coefficient_vapour = 0.; + cv.diffusion_coefficient_solute = 0.; + + cv.hCG = 0; + cv.hWG = 0; } } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h index 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/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..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, @@ -233,9 +231,9 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, // of the mass balance on the diagonal of the local element // matrix to be zero). The value is simply made up, seems // reasonable. - cv.xnWG = + double const xnWG = std::clamp(vapour_pressure_data.pWGR / pGR, xnWG_min, 1. - xnWG_min); - mass_mole_fractions_data.xnCG = 1. - cv.xnWG; + mass_mole_fractions_data.xnCG = 1. - xnWG; // gas phase molar fraction derivatives auto const dxnWG_dpGR = -vapour_pressure_data.pWGR / pGR / pGR; @@ -243,12 +241,12 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, auto const dxnWG_dT = dpWGR_dT / pGR; // molar mass of the gas phase as a mixture of 'air' and vapour - auto const MG = mass_mole_fractions_data.xnCG * M_C + cv.xnWG * M_W; + auto const MG = mass_mole_fractions_data.xnCG * M_C + xnWG * M_W; variables.molar_mass = MG; // gas phase mass fractions - cv.xmWG = cv.xnWG * M_W / MG; - mass_mole_fractions_data.xmCG = 1. - cv.xmWG; + double const xmWG = xnWG * M_W / MG; + mass_mole_fractions_data.xmCG = 1. - xmWG; auto const dxn_dxm_conversion = M_W * M_C / MG / MG; // gas phase mass fraction derivatives @@ -302,7 +300,7 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, // model is still consistent. constituent_density_data.rho_C_GR = mass_mole_fractions_data.xmCG * fluid_density_data.rho_GR; - constituent_density_data.rho_W_GR = cv.xmWG * fluid_density_data.rho_GR; + constituent_density_data.rho_W_GR = xmWG * fluid_density_data.rho_GR; // 'Air'-component partial density derivatives cv.drho_C_GR_dp_GR = mass_mole_fractions_data.xmCG * cv.drho_GR_dp_GR - @@ -314,11 +312,11 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, // Vapour-component partial density derivatives cv.drho_W_GR_dp_GR = - cv.xmWG * cv.drho_GR_dp_GR + cv.dxmWG_dpGR * fluid_density_data.rho_GR; - cv.drho_W_GR_dp_cap = cv.xmWG * cv.drho_GR_dp_cap + - cv.dxmWG_dpCap * fluid_density_data.rho_GR; + xmWG * cv.drho_GR_dp_GR + cv.dxmWG_dpGR * fluid_density_data.rho_GR; + cv.drho_W_GR_dp_cap = + xmWG * cv.drho_GR_dp_cap + cv.dxmWG_dpCap * fluid_density_data.rho_GR; cv.drho_W_GR_dT = - cv.xmWG * cv.drho_GR_dT + cv.dxmWG_dT * fluid_density_data.rho_GR; + xmWG * cv.drho_GR_dT + cv.dxmWG_dT * fluid_density_data.rho_GR; // specific heat capacities of dry air and vapour auto const cpCG = @@ -335,11 +333,13 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, cv.hWG = cpWG * T + dh_evap; // specific enthalpy of gas phase - enthalpy_data.h_G = - mass_mole_fractions_data.xmCG * cv.hCG + cv.xmWG * cv.hWG; + enthalpy_data.h_G = mass_mole_fractions_data.xmCG * cv.hCG + xmWG * cv.hWG; + cv.dh_G_dT = 0; // specific inner energies of gas phase cv.uG = enthalpy_data.h_G - pGR / fluid_density_data.rho_GR; + cv.du_G_dT = 0; + cv.du_G_dp_GR = 0; // diffusion auto const tortuosity = @@ -355,11 +355,6 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, variables.molar_fraction = mass_mole_fractions_data.xnCG; - // gas phase viscosity - viscosity_data.mu_GR = - gas_phase.property(MaterialPropertyLib::PropertyType::viscosity) - .template value(variables, x_t.x, x_t.t, x_t.dt); - // Dissolution part -- Liquid phase properties // ------------------------------------------- @@ -465,7 +460,7 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, // liquid phase density derivatives auto const drhoLR_dpGR = rho_ref_betaP + rho_ref_betaC * dcCL_dpGR; auto const drhoLR_dpCap = -rho_ref_betaP; - auto const drhoLR_dT = rho_ref_betaT + rho_ref_betaC * dcCL_dT; + cv.drho_LR_dT = rho_ref_betaT + rho_ref_betaC * dcCL_dT; // solvent partial density derivatives auto const drhoWLR_dpGR = rho_ref_betaP; @@ -480,7 +475,7 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, 1. / fluid_density_data.rho_LR * (drhoWLR_dpCap - mass_mole_fractions_data.xmWL * drhoLR_dpCap); cv.dxmWL_dT = 1. / fluid_density_data.rho_LR * - (drhoWLR_dT - mass_mole_fractions_data.xmWL * drhoLR_dT); + (drhoWLR_dT - mass_mole_fractions_data.xmWL * cv.drho_LR_dT); // liquid phase molar fractions and derivatives mass_mole_fractions_data.xnWL = @@ -508,12 +503,14 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, .template value(variables, x_t.x, x_t.t, x_t.dt); // specific enthalpy of liquid phase and its components - cv.hCL = cpCL * T + dh_sol; - cv.hWL = cpWL * T; - enthalpy_data.h_L = xmCL * cv.hCL + mass_mole_fractions_data.xmWL * cv.hWL; + double const hCL = cpCL * T + dh_sol; + double const hWL = cpWL * T; + enthalpy_data.h_L = xmCL * hCL + mass_mole_fractions_data.xmWL * hWL; + cv.dh_L_dT = 0; // specific inner energies of liquid phase cv.uL = enthalpy_data.h_L; + cv.du_L_dT = 0; // diffusion auto const D_C_L_m = @@ -522,10 +519,11 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, cv.diffusion_coefficient_solute = tortuosity * D_C_L_m; // Note here that D_C_L = D_W_L. - // liquid phase viscosity - viscosity_data.mu_LR = - liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity) - .template value(variables, x_t.x, x_t.t, x_t.dt); + // Some default initializations. + cv.drho_LR_dp_LR = 0; + cv.drho_W_LR_dp_GR = 0.; + cv.drho_W_LR_dT = 0.; + cv.drho_W_LR_dp_LR = 0.; } } // namespace ConstitutiveRelations diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h index 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/PhaseTransitionData.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h index 2500aa882ae..a0042c6afbb 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h @@ -17,68 +17,68 @@ namespace ConstitutiveRelations { struct PhaseTransitionData { - double drho_GR_dp_GR = 0.; - double drho_GR_dp_cap = 0.; - double drho_GR_dT = 0.; + double drho_GR_dp_GR = nan; + double drho_GR_dp_cap = nan; + double drho_GR_dT = nan; - double drho_C_GR_dp_GR = 0.; - double drho_C_GR_dp_cap = 0.; - double drho_C_GR_dT = 0.; + double drho_C_GR_dp_GR = nan; + double drho_C_GR_dp_cap = nan; + double drho_C_GR_dT = nan; - double drho_W_GR_dp_GR = 0.; - double drho_W_GR_dp_cap = 0.; - double drho_W_GR_dT = 0.; + double drho_W_GR_dp_GR = nan; + double drho_W_GR_dp_cap = nan; + double drho_W_GR_dT = nan; - double drho_LR_dp_GR = 0.; - double drho_LR_dp_cap = 0.; - double drho_LR_dT = 0.; - double drho_LR_dp_LR = 0.; + double drho_LR_dT = nan; + double drho_LR_dp_LR = nan; - double drho_C_LR_dp_GR = 0.; - double drho_C_LR_dp_cap = 0.; - double drho_C_LR_dT = 0.; - double drho_C_LR_dp_LR = 0.; + // TODO (naumov) These three are zero in both models but used in the + // assembly. Remove them and simplify assembly or correct the expressions in + // the phase transition models? + static constexpr double drho_C_LR_dp_GR = 0.; + static constexpr double drho_C_LR_dT = 0.; + static constexpr double drho_C_LR_dp_LR = 0.; - double drho_W_LR_dp_GR = 0.; - double drho_W_LR_dp_cap = 0.; - double drho_W_LR_dT = 0.; - double drho_W_LR_dp_LR = 0.; - - // constituent mass and molar fractions - double xnWG = 0.; - double xmWG = 0.; + double drho_W_LR_dp_GR = nan; + double drho_W_LR_dT = nan; + double drho_W_LR_dp_LR = nan; // mass fraction derivatives - double dxmWG_dpGR = 0.; - double dxmWG_dpCap = 0.; - double dxmWG_dT = 0.; + double dxmWG_dpGR = nan; + double dxmWG_dpCap = nan; + double dxmWG_dT = nan; - double dxmWL_dpGR = 0.; - double dxmWL_dpCap = 0.; - double dxmWL_dpLR = 0.; - double dxmWL_dT = 0.; + double dxmWL_dpGR = nan; + double dxmWL_dpCap = nan; + double dxmWL_dT = nan; + // TODO (naumov) This is zero in both models but used in the assembly. + // Remove it and simplify assembly or correct the expressions in the phase + // transition models? + static constexpr double dxmWL_dpLR = 0.; - double diffusion_coefficient_vapour = 0.; - double diffusion_coefficient_solute = 0.; + double diffusion_coefficient_vapour = nan; + double diffusion_coefficient_solute = nan; // specific enthalpies - double hCG = 0; - double hWG = 0; - double hCL = 0; - double hWL = 0; + double hCG = nan; + double hWG = nan; - double dh_G_dT = 0; - double dh_L_dT = 0; + double dh_G_dT = nan; + double dh_L_dT = nan; // specific inner energies - double uG = 0; - double uL = 0; - - double du_G_dT = 0; - double du_L_dT = 0; - double du_G_dp_GR = 0; - double du_L_dp_GR = 0; - double du_L_dp_cap = 0; + double uG = nan; + double uL = nan; + + double du_G_dT = nan; + double du_L_dT = nan; + double du_G_dp_GR = nan; + + // TODO (naumov) These two are zero in both models but used in the assembly. + // Remove them and simplify assembly or correct the expressions in the phase + // transition models? + static constexpr double du_L_dp_GR = 0; + static constexpr double du_L_dp_cap = 0; }; } // namespace ConstitutiveRelations diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h index 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/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/ConstitutiveRelations/SolidThermalExpansion.h b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h index 9f225f30f7d..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; /// Isotropic solid phase volumetric thermal expansion - /// coefficient. + double beta_T_SR = nan; /// Solid phase volumetric thermal expansion + /// coefficient. double thermal_volume_strain = nan; }; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp new file mode 100644 index 00000000000..abd9e0cc9f3 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp @@ -0,0 +1,112 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + */ + +#include "ThermalConductivity.h" + +#include "MaterialLib/MPL/Utils/FormEigenTensor.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +template +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; + + // Assuming dS_L/dT = 0, then: + // dphi_G_dT = -dS_L/dT * phi + (1 - S_L) * dphi_dT = (1 - S_L) * dphi_dT + // dphi_L_dT = dS_L/dT * phi + S_L * dphi_dT S_L * dphi_dT + // dphi_S_dT = -dphi_dT -dphi_dT + thermal_conductivity_data.dlambda_dT = + (1 - S_L_data.S_L) * porosity_d_data.dphi_dT * lambdaGR + + phi_G * dlambda_GR_dT + + S_L_data.S_L * porosity_d_data.dphi_dT * lambdaLR + + +phi_L * dlambda_LR_dT - porosity_d_data.dphi_dT * lambdaSR + + phi_S * dlambda_SR_dT; +} +template struct ThermalConductivityModel<2>; +template struct ThermalConductivityModel<3>; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h new file mode 100644 index 00000000000..fbf34021e55 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h @@ -0,0 +1,47 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + */ + +#pragma once + +#include "Base.h" +#include "Porosity.h" +#include "Saturation.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ + +template +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; +}; + +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/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/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h index ec9de5fdb91..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 = @@ -36,32 +29,7 @@ struct IntegrationPointData final typename ShapeMatricesTypePressure::NodalRowVectorType N_p; typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p; - // phase enthalpies - double rho_G_h_G = std::numeric_limits::quiet_NaN(); - 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(); - - GlobalDimMatrixType lambda; - GlobalDimVectorType d_CG; - GlobalDimVectorType d_WG; - GlobalDimVectorType d_CL; - GlobalDimVectorType d_WL; - - GlobalDimVectorType w_GS; - GlobalDimVectorType w_LS; - double integration_weight = std::numeric_limits::quiet_NaN(); - - void pushBackState() { rho_u_eff_prev = rho_u_eff; } - - EIGEN_MAKE_ALIGNED_OPERATOR_NEW; }; } // namespace TH2M diff --git a/ProcessLib/TH2M/LocalAssemblerInterface.h b/ProcessLib/TH2M/LocalAssemblerInterface.h index 9e8d86a0863..30e624d0b54 100644 --- a/ProcessLib/TH2M/LocalAssemblerInterface.h +++ b/ProcessLib/TH2M/LocalAssemblerInterface.h @@ -66,48 +66,6 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface, std::string_view name, double const* values, int const integration_order) = 0; - virtual std::vector 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; - - 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; - - 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 9203a3d5326..2f14daaaeeb 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -115,9 +115,6 @@ 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 = @@ -162,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; @@ -227,11 +223,15 @@ TH2MLocalAssemblerS_L; - - vars.porosity = ip_out.porosity_data.phi; + models.thermal_conductivity_model.eval( + {pos, t, dt}, media_data, T_data, ip_out.porosity_data, + ip_cv.porosity_d_data, current_state.S_L_data, ip_cv.dS_L_dp_cap, + ip_cv.thermal_conductivity_data); auto const& c = ip_cv.phase_transition_data; @@ -269,29 +262,20 @@ TH2MLocalAssembler( - medium - .property( - MaterialPropertyLib::PropertyType::thermal_conductivity) - .value(vars, pos, t, dt)); - - auto const cpS = - solid_phase.property(MPL::PropertyType::specific_heat_capacity) - .template value(vars, pos, t, dt); - ip_data.h_S = cpS * T; - auto const u_S = ip_data.h_S; - - ip_data.rho_u_eff = phi_G * ip_out.fluid_density_data.rho_GR * c.uG + - phi_L * ip_out.fluid_density_data.rho_LR * c.uL + - phi_S * ip_out.solid_density_data.rho_SR * u_S; - - ip_data.rho_G_h_G = - phi_G * ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G; - ip_data.rho_L_h_L = - phi_L * ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L; - ip_data.rho_S_h_S = - phi_S * ip_out.solid_density_data.rho_SR * ip_data.h_S; + ip_out.enthalpy_data.h_S = ip_cv.solid_heat_capacity_data() * T; + auto const u_S = ip_out.enthalpy_data.h_S; + + current_state.internal_energy_data() = + phi_G * ip_out.fluid_density_data.rho_GR * c.uG + + phi_L * ip_out.fluid_density_data.rho_LR * c.uL + + phi_S * ip_out.solid_density_data.rho_SR * u_S; + + ip_cv.effective_volumetric_enthalpy_data.rho_h_eff = + phi_G * ip_out.fluid_density_data.rho_GR * + ip_out.enthalpy_data.h_G + + phi_L * ip_out.fluid_density_data.rho_LR * + ip_out.enthalpy_data.h_L + + phi_S * ip_out.solid_density_data.rho_SR * ip_out.enthalpy_data.h_S; // for variable output auto const xmCL = 1. - ip_out.mass_mole_fractions_data.xmWL; @@ -310,45 +294,44 @@ TH2MLocalAssembler(vars, MPL::Variable::temperature, pos, - t, dt); - - ip_cv.drho_u_eff_dT = + ip_cv.effective_volumetric_internal_energy_d_data.drho_u_eff_dT = phi_G * c.drho_GR_dT * c.uG + phi_G * ip_out.fluid_density_data.rho_GR * c.du_G_dT + - phi_L * drho_LR_dT * c.uL + + phi_L * c.drho_LR_dT * c.uL + phi_L * ip_out.fluid_density_data.rho_LR * c.du_L_dT + phi_S * ip_cv.solid_density_d_data.drho_SR_dT * u_S + - phi_S * ip_out.solid_density_data.rho_SR * cpS - + phi_S * ip_out.solid_density_data.rho_SR * + ip_cv.solid_heat_capacity_data() - ip_cv.porosity_d_data.dphi_dT * ip_out.solid_density_data.rho_SR * u_S; @@ -363,58 +346,6 @@ TH2MLocalAssembler( - 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.dlambda_dp_cap = - dphi_G_dp_cap * lambdaGR + dphi_L_dp_cap * lambdaLR; - - ip_cv.dlambda_dT = phi_G * dlambda_GR_dT + phi_L * dlambda_LR_dT + - phi_S * dlambda_SR_dT - - ip_cv.porosity_d_data.dphi_dT * lambdaSR; - // From p_LR = p_GR - p_cap it follows for // drho_LR/dp_GR = drho_LR/dp_LR * dp_LR/dp_GR // = drho_LR/dp_LR * (dp_GR/dp_GR - dp_cap/dp_GR) @@ -423,14 +354,14 @@ TH2MLocalAssemblermaterial_states_[ip].pushBackState(); this->prev_states_[ip] = this->current_states_[ip]; } @@ -1199,10 +1137,6 @@ void TH2MLocalAssembler< auto const rho_W_LR_dot = (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt; - auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S; - - auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt; - GlobalDimMatrixType const k_over_mu_G = ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G / ip_cv.viscosity_data.mu_GR; @@ -1388,27 +1322,38 @@ void TH2MLocalAssembler< // - temperature equation // --------------------------------------------------------------------- - MTu.noalias() += NTT * rho_h_eff * mT * Bu * w; + MTu.noalias() += NTT * + ip_cv.effective_volumetric_enthalpy_data.rho_h_eff * + mT * Bu * w; - KTT.noalias() += gradNTT * ip.lambda * gradNT * w; + KTT.noalias() += + gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w; + auto const rho_u_eff_dot = (current_state.internal_energy_data() - + **prev_state.internal_energy_data) / + dt; fT.noalias() -= NTT * rho_u_eff_dot * w; fT.noalias() += gradNTT * - (rhoGR * ip_out.enthalpy_data.h_G * ip.w_GS + - rhoLR * ip_out.enthalpy_data.h_L * ip.w_LS) * + (rhoGR * ip_out.enthalpy_data.h_G * + ip_out.darcy_velocity_data.w_GS + + rhoLR * ip_out.enthalpy_data.h_L * + ip_out.darcy_velocity_data.w_LS) * w; fT.noalias() += gradNTT * (current_state.constituent_density_data.rho_C_GR * - ip_cv.phase_transition_data.hCG * ip.d_CG + + ip_cv.phase_transition_data.hCG * + ip_out.diffusion_velocity_data.d_CG + current_state.constituent_density_data.rho_W_GR * - ip_cv.phase_transition_data.hWG * ip.d_WG) * + ip_cv.phase_transition_data.hWG * + ip_out.diffusion_velocity_data.d_WG) * w; - fT.noalias() += - NTT * (rhoGR * ip.w_GS.transpose() + rhoLR * ip.w_LS.transpose()) * - b * w; + fT.noalias() += NTT * + (rhoGR * ip_out.darcy_velocity_data.w_GS.transpose() + + rhoLR * ip_out.darcy_velocity_data.w_LS.transpose()) * + b * w; // --------------------------------------------------------------------- // - displacement equation @@ -1690,10 +1635,6 @@ void TH2MLocalAssembler(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.lambda * gradNT * w; + KTT.noalias() += + gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w; // d KTT/dp_GR * T // TODO (naumov) always zero if lambda_xR have no derivatives wrt. p_GR. @@ -2071,40 +2021,57 @@ void TH2MLocalAssembler(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 + auto const rho_u_eff_dot = (current_state.internal_energy_data() - + **prev_state.internal_energy_data) / + dt; fT.noalias() -= NTT * rho_u_eff_dot * w; // dfT_1/dp_GR local_Jac .template block(temperature_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 * - (rhoGR * ip_out.enthalpy_data.h_G * ip.w_GS + - rhoLR * ip_out.enthalpy_data.h_L * ip.w_LS) * + (rhoGR * ip_out.enthalpy_data.h_G * + ip_out.darcy_velocity_data.w_GS + + rhoLR * ip_out.enthalpy_data.h_L * + ip_out.darcy_velocity_data.w_LS) * w; // dfT_2/dp_GR @@ -2134,15 +2101,18 @@ void TH2MLocalAssembler -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 -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 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(); } @@ -190,37 +188,6 @@ class TH2MLocalAssembler : public LocalAssemblerInterface return Eigen::Map(N_u.data(), N_u.size()); } - 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::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>, std::vector< @@ -231,28 +198,12 @@ 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; - using IpData = - IntegrationPointData; - std::vector> _ip_data; + using IpData = IntegrationPointData; + std::vector _ip_data; SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType> diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp index 34d2e0d1db4..b3289ca9d14 100644 --- a/ProcessLib/TH2M/TH2MProcess.cpp +++ b/ProcessLib/TH2M/TH2MProcess.cpp @@ -163,33 +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); - add_secondary_variable( - "diffusion_velocity_vapour_gas", mesh.getDimension(), - &LocalAssemblerInterface< - DisplacementDim>::getIntPtDiffusionVelocityVapourGas); - add_secondary_variable( - "diffusion_velocity_gas_gas", mesh.getDimension(), - &LocalAssemblerInterface< - DisplacementDim>::getIntPtDiffusionVelocityGasGas); - add_secondary_variable( - "diffusion_velocity_solute_liquid", mesh.getDimension(), - &LocalAssemblerInterface< - DisplacementDim>::getIntPtDiffusionVelocitySoluteLiquid); - add_secondary_variable( - "diffusion_velocity_liquid_liquid", mesh.getDimension(), - &LocalAssemblerInterface< - DisplacementDim>::getIntPtDiffusionVelocityLiquidLiquid); - - add_secondary_variable( - "enthalpy_solid", 1, - &LocalAssemblerInterface::getIntPtEnthalpySolid); - // // enable output of internal variables defined by material models // 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 diff --git a/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp b/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp index 52b49cd7a02..bc74bda9cf0 100644 --- a/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp +++ b/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp @@ -34,9 +34,6 @@ TEST(ProcessLib, TH2MNoPhaseTransition) auto const diffusion_coefficient_vapour = 0.0; - auto const viscosity_air = 1.e-5; - auto const viscosity_water = 1.e-3; - m << "\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,17 +135,14 @@ TEST(ProcessLib, TH2MNoPhaseTransition) double const hL = specific_heat_capacity_water * T; double const uG = hG - pGR / density_air; double const uL = hL; - double const muGR = viscosity_air; - double const muLR = viscosity_water; ASSERT_NEAR(density_air, fluid_density.rho_GR, 1.0e-10); ASSERT_NEAR(density_water, fluid_density.rho_LR, 1.0e-10); ASSERT_NEAR(rhoCGR, constituent_density.rho_C_GR, 1.0e-10); ASSERT_NEAR(rhoWGR, constituent_density.rho_W_GR, 1.0e-10); ASSERT_NEAR(rhoCLR, constituent_density.rho_C_LR, 1.0e-10); - ASSERT_NEAR(xmCG, 1. - cv.xmWG, 1.e-10); ASSERT_NEAR(xmCG, mass_mole_fractions.xmCG, 1.e-10); - ASSERT_NEAR(xmWG, cv.xmWG, 1.e-10); + ASSERT_NEAR(xmWG, 1 - mass_mole_fractions.xmCG, 1.e-10); ASSERT_NEAR(dxmWG_dpGR, cv.dxmWG_dpGR, 1.0e-10); ASSERT_NEAR(dxmCG_dpGR, -cv.dxmWG_dpGR, 1.0e-10); ASSERT_NEAR(dxmWG_dT, cv.dxmWG_dT, 1.0e-10); @@ -164,6 +155,4 @@ TEST(ProcessLib, TH2MNoPhaseTransition) ASSERT_NEAR(uL, cv.uL, 1.0e-10); ASSERT_NEAR(diffusion_coefficient_vapour, cv.diffusion_coefficient_vapour, 1.0e-10); - ASSERT_NEAR(muGR, viscosity.mu_GR, 1.0e-10); - ASSERT_NEAR(muLR, viscosity.mu_LR, 1.0e-10); } diff --git a/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp b/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp index 7fbe098c97d..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,11 +241,10 @@ TEST(ProcessLib, TH2MPhaseTransition) rhoWLR_model.eval(x_t, media_data, CR::GasPressureData{pGR + eps_pGR}, p_cap_data, T_data, rhoWLR); ptm->eval(x_t, media_data, CR::GasPressureData{pGR + eps_pGR}, - p_cap_data, T_data, rhoWLR, viscosity, enthalpy, - mass_mole_fractions, fluid_density, vapour_pressure, - constituent_density, cv); + p_cap_data, T_data, rhoWLR, enthalpy, mass_mole_fractions, + fluid_density, vapour_pressure, constituent_density, cv); - auto xmWG_plus = cv.xmWG; + auto xmWG_plus = 1 - mass_mole_fractions.xmCG; auto rhoGR_plus = fluid_density.rho_GR; auto rhoCGR_plus = constituent_density.rho_C_GR; auto rhoWGR_plus = constituent_density.rho_W_GR; @@ -258,11 +252,10 @@ TEST(ProcessLib, TH2MPhaseTransition) rhoWLR_model.eval(x_t, media_data, CR::GasPressureData{pGR - eps_pGR}, p_cap_data, T_data, rhoWLR); ptm->eval(x_t, media_data, CR::GasPressureData{pGR - eps_pGR}, - p_cap_data, T_data, rhoWLR, viscosity, enthalpy, - mass_mole_fractions, fluid_density, vapour_pressure, - constituent_density, cv); + p_cap_data, T_data, rhoWLR, enthalpy, mass_mole_fractions, + fluid_density, vapour_pressure, constituent_density, cv); - auto xmWG_minus = cv.xmWG; + auto xmWG_minus = 1 - mass_mole_fractions.xmCG; auto rhoGR_minus = fluid_density.rho_GR; auto rhoCGR_minus = constituent_density.rho_C_GR; auto rhoWGR_minus = constituent_density.rho_W_GR; @@ -271,8 +264,8 @@ TEST(ProcessLib, TH2MPhaseTransition) rhoWLR_model.eval(x_t, media_data, CR::GasPressureData{pGR}, p_cap_data, T_data, rhoWLR); ptm->eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR, - viscosity, enthalpy, mass_mole_fractions, fluid_density, - vapour_pressure, constituent_density, cv); + enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + constituent_density, cv); // Central difference derivatives auto const dxmWG_dpGR = (xmWG_plus - xmWG_minus) / (2. * eps_pGR); @@ -297,10 +290,10 @@ TEST(ProcessLib, TH2MPhaseTransition) rhoWLR); ptm->eval(x_t, media_data, p_GR_data, CR::CapillaryPressureData{pCap + eps_pCap}, T_data, rhoWLR, - viscosity, enthalpy, mass_mole_fractions, fluid_density, - vapour_pressure, constituent_density, cv); + enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + constituent_density, cv); - xmWG_plus = cv.xmWG; + xmWG_plus = 1 - mass_mole_fractions.xmCG; rhoGR_plus = fluid_density.rho_GR; rhoCGR_plus = constituent_density.rho_C_GR; rhoWGR_plus = constituent_density.rho_W_GR; @@ -310,10 +303,10 @@ TEST(ProcessLib, TH2MPhaseTransition) rhoWLR); ptm->eval(x_t, media_data, p_GR_data, CR::CapillaryPressureData{pCap - eps_pCap}, T_data, rhoWLR, - viscosity, enthalpy, mass_mole_fractions, fluid_density, - vapour_pressure, constituent_density, cv); + enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + constituent_density, cv); - xmWG_minus = cv.xmWG; + xmWG_minus = 1 - mass_mole_fractions.xmCG; rhoGR_minus = fluid_density.rho_GR; rhoCGR_minus = constituent_density.rho_C_GR; rhoWGR_minus = constituent_density.rho_W_GR; @@ -322,8 +315,8 @@ TEST(ProcessLib, TH2MPhaseTransition) rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR); ptm->eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR, - viscosity, enthalpy, mass_mole_fractions, fluid_density, - vapour_pressure, constituent_density, cv); + enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + constituent_density, cv); // Central difference derivatives auto const dxmWG_dpCap = (xmWG_plus - xmWG_minus) / (2. * eps_pCap); @@ -348,11 +341,11 @@ TEST(ProcessLib, TH2MPhaseTransition) rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, CR::TemperatureData{T + eps_T, T}, rhoWLR); ptm->eval(x_t, media_data, p_GR_data, p_cap_data, - CR::TemperatureData{T + eps_T, T}, rhoWLR, viscosity, - enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + CR::TemperatureData{T + eps_T, T}, rhoWLR, enthalpy, + mass_mole_fractions, fluid_density, vapour_pressure, constituent_density, cv); - xmWG_plus = cv.xmWG; + xmWG_plus = 1 - mass_mole_fractions.xmCG; rhoGR_plus = fluid_density.rho_GR; rhoCGR_plus = constituent_density.rho_C_GR; rhoWGR_plus = constituent_density.rho_W_GR; @@ -360,11 +353,11 @@ TEST(ProcessLib, TH2MPhaseTransition) rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, CR::TemperatureData{T - eps_T, T}, rhoWLR); ptm->eval(x_t, media_data, p_GR_data, p_cap_data, - CR::TemperatureData{T - eps_T, T}, rhoWLR, viscosity, - enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + CR::TemperatureData{T - eps_T, T}, rhoWLR, enthalpy, + mass_mole_fractions, fluid_density, vapour_pressure, constituent_density, cv); - xmWG_minus = cv.xmWG; + xmWG_minus = 1 - mass_mole_fractions.xmCG; rhoGR_minus = fluid_density.rho_GR; rhoCGR_minus = constituent_density.rho_C_GR; rhoWGR_minus = constituent_density.rho_W_GR; @@ -373,8 +366,8 @@ TEST(ProcessLib, TH2MPhaseTransition) rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR); ptm->eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR, - viscosity, enthalpy, mass_mole_fractions, fluid_density, - vapour_pressure, constituent_density, cv); + enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + constituent_density, cv); // Central difference derivatives auto const dxmWG_dT = (xmWG_plus - xmWG_minus) / (2. * eps_T); @@ -397,18 +390,19 @@ TEST(ProcessLib, TH2MPhaseTransition) // phase pressure auto const refence_xnWG = std::clamp(vapour_pressure.pWGR / pGR, 0., 1.); - ASSERT_NEAR(refence_xnWG, cv.xnWG, 1.e-10); + auto const xnWG = 1. - mass_mole_fractions.xnCG; + ASSERT_NEAR(refence_xnWG, xnWG, 1.e-10); // The quotient of constituent partial densities and phase densities // must be equal to the mass fraction of those constituents in both // phases. ASSERT_NEAR(constituent_density.rho_W_GR / fluid_density.rho_GR, - cv.xmWG, 1.e-10); + 1 - mass_mole_fractions.xmCG, 1.e-10); ASSERT_NEAR(rhoWLR() / fluid_density.rho_LR, mass_mole_fractions.xmWL, 1.e-10); ASSERT_NEAR(constituent_density.rho_C_GR / fluid_density.rho_GR, - 1. - cv.xmWG, 1.e-10); + mass_mole_fractions.xmCG, 1.e-10); ASSERT_NEAR(constituent_density.rho_C_LR / fluid_density.rho_LR, 1. - mass_mole_fractions.xmWL, 1.e-10); @@ -424,8 +418,8 @@ TEST(ProcessLib, TH2MPhaseTransition) // Gas density (ideal gas in this test): constexpr double R = MaterialLib::PhysicalConstant::IdealGasConstant; - auto const xnCG = 1. - cv.xnWG; - auto const MG = cv.xnWG * molar_mass_water + xnCG * molar_mass_air; + auto const MG = + xnWG * molar_mass_water + mass_mole_fractions.xnCG * molar_mass_air; auto const rhoGR = pGR * MG / R / T; ASSERT_NEAR(rhoGR, fluid_density.rho_GR, 1.e-10); @@ -449,9 +443,6 @@ TEST(ProcessLib, TH2MPhaseTransition) diffusion_coefficient_vapour; ASSERT_NEAR(diffusion_gas_phase, cv.diffusion_coefficient_vapour, 1.0e-10); - - ASSERT_NEAR(viscosity_air, viscosity.mu_GR, 1.0e-10); - ASSERT_NEAR(viscosity_water, viscosity.mu_LR, 1.0e-10); } } @@ -489,7 +480,6 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho) CR::PureLiquidDensityData rhoWLR; CR::PureLiquidDensityModel rhoWLR_model; - CR::ViscosityData viscosity; CR::EnthalpyData enthalpy; CR::MassMoleFractionsData mass_mole_fractions; CR::FluidDensityData fluid_density; @@ -513,22 +503,20 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho) rhoWLR_model.eval(x_t, media_data, CR::GasPressureData{pGR + eps_pGR}, p_cap_data, T_data, rhoWLR); ptm->eval(x_t, media_data, CR::GasPressureData{pGR + eps_pGR}, - p_cap_data, T_data, rhoWLR, viscosity, enthalpy, - mass_mole_fractions, fluid_density, vapour_pressure, - constituent_density, cv); + p_cap_data, T_data, rhoWLR, enthalpy, mass_mole_fractions, + fluid_density, vapour_pressure, constituent_density, cv); - auto xmWG_plus = cv.xmWG; + auto xmWG_plus = 1 - mass_mole_fractions.xmCG; auto rhoCGR_plus = constituent_density.rho_C_GR; auto rhoWGR_plus = constituent_density.rho_W_GR; rhoWLR_model.eval(x_t, media_data, CR::GasPressureData{pGR - eps_pGR}, p_cap_data, T_data, rhoWLR); ptm->eval(x_t, media_data, CR::GasPressureData{pGR - eps_pGR}, - p_cap_data, T_data, rhoWLR, viscosity, enthalpy, - mass_mole_fractions, fluid_density, vapour_pressure, - constituent_density, cv); + p_cap_data, T_data, rhoWLR, enthalpy, mass_mole_fractions, + fluid_density, vapour_pressure, constituent_density, cv); - auto xmWG_minus = cv.xmWG; + auto xmWG_minus = 1 - mass_mole_fractions.xmCG; auto rhoCGR_minus = constituent_density.rho_C_GR; auto rhoWGR_minus = constituent_density.rho_W_GR; @@ -536,8 +524,8 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho) rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR); ptm->eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR, - viscosity, enthalpy, mass_mole_fractions, fluid_density, - vapour_pressure, constituent_density, cv); + enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + constituent_density, cv); // Central difference derivatives auto const dxmWG_dpGR = (xmWG_plus - xmWG_minus) / (2. * eps_pGR); @@ -561,10 +549,10 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho) rhoWLR); ptm->eval(x_t, media_data, p_GR_data, CR::CapillaryPressureData{pCap + eps_pCap}, T_data, rhoWLR, - viscosity, enthalpy, mass_mole_fractions, fluid_density, - vapour_pressure, constituent_density, cv); + enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + constituent_density, cv); - xmWG_plus = cv.xmWG; + xmWG_plus = 1 - mass_mole_fractions.xmCG; rhoCGR_plus = constituent_density.rho_C_GR; rhoWGR_plus = constituent_density.rho_W_GR; @@ -573,10 +561,10 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho) rhoWLR); ptm->eval(x_t, media_data, p_GR_data, CR::CapillaryPressureData{pCap - eps_pCap}, T_data, rhoWLR, - viscosity, enthalpy, mass_mole_fractions, fluid_density, - vapour_pressure, constituent_density, cv); + enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + constituent_density, cv); - xmWG_minus = cv.xmWG; + xmWG_minus = 1 - mass_mole_fractions.xmCG; rhoCGR_minus = constituent_density.rho_C_GR; rhoWGR_minus = constituent_density.rho_W_GR; @@ -584,8 +572,8 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho) rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR); ptm->eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR, - viscosity, enthalpy, mass_mole_fractions, fluid_density, - vapour_pressure, constituent_density, cv); + enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + constituent_density, cv); // Central difference derivatives auto const dxmWG_dpCap = (xmWG_plus - xmWG_minus) / (2. * eps_pCap); @@ -609,22 +597,22 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho) rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, CR::TemperatureData{T + eps_T, T}, rhoWLR); ptm->eval(x_t, media_data, p_GR_data, p_cap_data, - CR::TemperatureData{T + eps_T, T}, rhoWLR, viscosity, - enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + CR::TemperatureData{T + eps_T, T}, rhoWLR, enthalpy, + mass_mole_fractions, fluid_density, vapour_pressure, constituent_density, cv); - xmWG_plus = cv.xmWG; + xmWG_plus = 1 - mass_mole_fractions.xmCG; rhoCGR_plus = constituent_density.rho_C_GR; rhoWGR_plus = constituent_density.rho_W_GR; rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, CR::TemperatureData{T - eps_T, T}, rhoWLR); ptm->eval(x_t, media_data, p_GR_data, p_cap_data, - CR::TemperatureData{T - eps_T, T}, rhoWLR, viscosity, - enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + CR::TemperatureData{T - eps_T, T}, rhoWLR, enthalpy, + mass_mole_fractions, fluid_density, vapour_pressure, constituent_density, cv); - xmWG_minus = cv.xmWG; + xmWG_minus = 1 - mass_mole_fractions.xmCG; rhoCGR_minus = constituent_density.rho_C_GR; rhoWGR_minus = constituent_density.rho_W_GR; @@ -632,8 +620,8 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho) rhoWLR_model.eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR); ptm->eval(x_t, media_data, p_GR_data, p_cap_data, T_data, rhoWLR, - viscosity, enthalpy, mass_mole_fractions, fluid_density, - vapour_pressure, constituent_density, cv); + enthalpy, mass_mole_fractions, fluid_density, vapour_pressure, + constituent_density, cv); // Central difference derivatives auto const dxmWG_dT = (xmWG_plus - xmWG_minus) / (2. * eps_T); @@ -655,18 +643,19 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho) // phase pressure auto const refence_xnWG = std::clamp(vapour_pressure.pWGR / pGR, 0., 1.); - ASSERT_NEAR(refence_xnWG, cv.xnWG, 1.e-10); + auto const xnWG = 1. - mass_mole_fractions.xnCG; + ASSERT_NEAR(refence_xnWG, xnWG, 1.e-10); // The quotient of constituent partial densities and phase densities // must be equal to the mass fraction of those constituents in both // phases. ASSERT_NEAR(constituent_density.rho_W_GR / fluid_density.rho_GR, - cv.xmWG, 1.e-10); + 1 - mass_mole_fractions.xmCG, 1.e-10); ASSERT_NEAR(rhoWLR() / fluid_density.rho_LR, mass_mole_fractions.xmWL, 1.e-10); ASSERT_NEAR(constituent_density.rho_C_GR / fluid_density.rho_GR, - 1. - cv.xmWG, 1.e-10); + mass_mole_fractions.xmCG, 1.e-10); ASSERT_NEAR(constituent_density.rho_C_LR / fluid_density.rho_LR, 1. - mass_mole_fractions.xmWL, 1.e-10); @@ -704,8 +693,5 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho) diffusion_coefficient_vapour; ASSERT_NEAR(diffusion_gas_phase, cv.diffusion_coefficient_vapour, 1.0e-10); - - ASSERT_NEAR(viscosity_air, viscosity.mu_GR, 1.0e-10); - ASSERT_NEAR(viscosity_water, viscosity.mu_LR, 1.0e-10); } }