From 68ebbeccd1076db2bc17f2ce2e8576e74cfe2909 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Tue, 22 Nov 2022 18:32:30 +0100
Subject: [PATCH] [PL/THM] Implement freezing for temperature eq.
Using volume_fraction and frozen liquid properties.
Avoiding SpecificHeatCapacityWithLatentHeat; the approach
results in quantities not appearing in the PDE in particular
the quotient of effective volumetric specific heat capacity
and effective density (C_vol/rho_eff).
---
.../ThermoHydroMechanicsFEM-impl.h | 67 +++++++++++++++++--
1 file changed, 62 insertions(+), 5 deletions(-)
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
index 009cf7ca9b9..a91ac8f4b79 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
@@ -229,6 +229,9 @@ void ThermoHydroMechanicsLocalAssembler<
auto const& medium = _process_data.media_map->getMedium(_element.getID());
auto const& liquid_phase = medium->phase("AqueousLiquid");
auto const& solid_phase = medium->phase("Solid");
+ auto* const frozen_liquid_phase = medium->hasPhase("FrozenLiquid")
+ ? &medium->phase("FrozenLiquid")
+ : nullptr;
MaterialPropertyLib::VariableArray vars;
typename ShapeMatricesTypePressure::NodalVectorType node_flux_q;
@@ -484,14 +487,68 @@ void ThermoHydroMechanicsLocalAssembler<
std::max(max_velocity_magnitude, velocity.norm());
}
- auto const effective_volumetric_heat_capacity =
+ double const c_s =
+ solid_phase
+ .property(
+ MaterialPropertyLib::PropertyType::specific_heat_capacity)
+ .template value<double>(vars, x_position, t, dt);
+
+ // Non-const for modification by freezing terms.
+ double effective_volumetric_heat_capacity =
porosity * fluid_density * c_f +
- (1.0 - porosity) * solid_density *
- solid_phase
- .property(MaterialPropertyLib::PropertyType::
- specific_heat_capacity)
+ (1.0 - porosity) * solid_density * c_s;
+
+ if (frozen_liquid_phase)
+ {
+ double const phi_fr =
+ (*medium)[MaterialPropertyLib::PropertyType::volume_fraction]
.template value<double>(vars, x_position, t, dt);
+ auto const frozen_liquid_value =
+ [&](MaterialPropertyLib::PropertyType const p)
+ {
+ return (*frozen_liquid_phase)[p].template value<double>(
+ vars, x_position, t, dt);
+ };
+
+ double const rho_fr =
+ frozen_liquid_value(MaterialPropertyLib::PropertyType::density);
+
+ double const c_fr = frozen_liquid_value(
+ MaterialPropertyLib::PropertyType::specific_heat_capacity);
+
+ // specific latent heat of melting
+ double const l_fr = frozen_liquid_value(
+ MaterialPropertyLib::PropertyType::specific_latent_heat);
+
+ double const dphi_fr_dT =
+ (*medium)[MaterialPropertyLib::PropertyType::volume_fraction]
+ .template dValue<double>(
+ vars, MaterialPropertyLib::Variable::temperature,
+ x_position, t, dt);
+
+ effective_volumetric_heat_capacity +=
+ -phi_fr * fluid_density * c_f + phi_fr * rho_fr * c_fr -
+ l_fr * rho_fr * dphi_fr_dT;
+
+ // part of dMTT_dT derivative for freezing
+ double const d2phi_fr_dT2 =
+ (*medium)[MaterialPropertyLib::PropertyType::volume_fraction]
+ .template d2Value<double>(
+ vars, MaterialPropertyLib::Variable::temperature,
+ MaterialPropertyLib::Variable::temperature, x_position,
+ t, dt);
+
+ local_Jac
+ .template block<temperature_size, temperature_size>(
+ temperature_index, temperature_index)
+ .noalias() -=
+ N_T.transpose() *
+ ((rho_fr * c_fr - fluid_density * c_f) * dphi_fr_dT +
+ l_fr * rho_fr * d2phi_fr_dT2) *
+ N_T.dot(T_dot) * N_T * w;
+ }
+
MTT.noalias() +=
N_T.transpose() * effective_volumetric_heat_capacity * N_T * w;
--
GitLab