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.

Author: andrewdnolan

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