From a6041dbb8da42aa67f6c5ae5575866a24d9e35ac Mon Sep 17 00:00:00 2001 From: andrewdnolan Date: Tue, 3 May 2022 16:46:45 -0700 Subject: [PATCH] switched from Gilbert et al. 2014 to Zwinger et al. 2007 for thermal conducivity of surface layer to account for both temp and density effects. --- src/elmer_UDF/Thermodynamics.f90 | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/elmer_UDF/Thermodynamics.f90 b/src/elmer_UDF/Thermodynamics.f90 index 1460125..bd533ed 100644 --- a/src/elmer_UDF/Thermodynamics.f90 +++ b/src/elmer_UDF/Thermodynamics.f90 @@ -6,19 +6,19 @@ function Diffusivity(Model, Node, Temp) result(EnthalpyDiffusivity) type(Variable_t), pointer :: Density integer :: Node ! current node number - real(kind=dp) :: Temp, & ! [K] - T_ptr, & ! [K] Temp. (of water) at triple point - rho, & ! [kg m^-3] + real(kind=dp) :: Temp, & ! [K] + rho, & ! [kg m^-3] + rho_i, & ! [kg m^-3] Heat_Capacity, & ! [J kg^-1 K^-1 EnthalpyDiffusivity, & ! [kg m^-1 a^-1] HeatConductivity, & ! [W K^-1 m^-1] thermal conductivity - K_rho, & ! [W K^-1 m^-1] desnity dependence of "" - K_ice, & ! [W K^-1 m^-1] "" of ice - K_ptr, & ! [W K^-1 m^-1] "" at triple point of water - CapA, & ! [J kg-1 K-2] - CapB, & ! [J kg-1 K-1] - CondA, & ! [W m^5 K^-1 kg^-2] - CondB, & ! [W m^2 K^-1 kg^-1] + K_rho, & ! [W K^-1 m^-1] desnity dependence of "" + K_ice, & ! [W K^-1 m^-1] "" of ice + K_rho_i, & ! [W K^-1 m^-1] "" at rho of ice + CapA, & ! [J kg-1 K-2] + CapB, & ! [J kg-1 K-1] + CondA, & ! [W m^5 K^-1 kg^-2] + CondB, & ! [W m^2 K^-1 kg^-1] CondC ! [W kg-1 m-1] Density => VariableGet(Model % Variables, 'Densi') ! [kg m^-3] @@ -27,7 +27,7 @@ function Diffusivity(Model, Node, Temp) result(EnthalpyDiffusivity) CondA = GetConstReal(Model % Constants, "Enthalpy Heat Conductivity A") ! [W m^5 K^-1 kg^-2] CondB = GetConstReal(Model % Constants, "Enthalpy Heat Conductivity B") ! [W m^2 K^-1 kg^-1] CondC = GetConstReal(Model % Constants, "Enthalpy Heat Conductivity C") ! [W kg-1 m-1] - T_ptr = GetConstReal(Model % Constants, "T_triple") ! [K] + rho_i = GetConstReal(Model % Constants, "rho_i") ! [K] Temp = Temp + 273.15 ! [K] <-- [C] rho = Density%values(Density%perm(Node)) ! [kg m^-3] @@ -36,12 +36,12 @@ function Diffusivity(Model, Node, Temp) result(EnthalpyDiffusivity) Heat_Capacity = CapA * Temp + CapB ! Intermediate conductivity calcs. [W kg-1 m-1]] - K_rho = CondA * rho**2 - CondB * rho + CondC + K_rho = CondA * rho**2 - CondB * rho + CondC + K_rho_i = CondA * rho_i**2 - CondB * rho_i + CondC K_ice = 9.828*exp(-5.7e-3 *(Temp)) - K_ptr = 9.828*exp(-5.7e-3 *(T_ptr)) ! Conductivity [J a^-1 kg^-1 K^-1] <-- [W K^-1 m^-1] == [J s^-1 K^-1 m^-1] - HeatConductivity = (K_ice / K_ptr * K_rho) * 3600.0*24.0*365.25 + HeatConductivity = (K_rho / K_rho_i * K_ice) * 3600.0*24.0*365.25 ! Diffusivity [kg m^-1 a^-1] EnthalpyDiffusivity = HeatConductivity / Heat_Capacity