Merged PR

python/KaempferPlapp2009.py

@@ -35,15 +35,29 @@
 R_da = 286.9
 R_v = 461.5
 P_a = 1.01325e5
+rho_a = 1.341
+rho_i = 918.9
 K_fit = [-0.58653696e4, 0.2224103300e2, 0.13749042e-1, -0.34031775e-4, 0.26967687e-7, 0.6918651]
 
 # Eq. (2)
 P_vs = exp(K_fit[0]*T**(-1) + K_fit[1]*T**0 + K_fit[2]*T**1 + K_fit[3]*T**2 + K_fit[4]*T**3 + K_fit[5]*log(T))
 
-# THIS IS GIVING THE WRONG NUMBER!
 # Eq. (1)
 x_s = (R_da / R_v) * P_vs / (P_a - P_vs)
 
+# Eq. (3)
+rho_vs = rho_a * x_s
+
+# Concentration equilibrium
+u_eq = (rho_vs - rho_vs.subs(T, 263.15)) / rho_i
+
+
 print "T = 263.15"
 print "P_vs = ", P_vs.evalf(subs={T: 263.15})
 print "x_s = ", x_s.evalf(subs={T: 263.15})
+print "rho_vs(263.15) = ", rho_vs.evalf(subs={T: 263.15})
+print "rho_vs(268.15) = ", rho_vs.evalf(subs={T: 263.15})
+
+print ""
+print "u_eq(263.15) = ", u_eq.evalf(subs={T: 263.15})
+print "u_eq(268.15) = ", u_eq.evalf(subs={T: 268.15})

src/materials/PhaseFieldProperties.C

@@ -93,6 +93,6 @@ PhaseFieldProperties::computeQpProperties()
   // x_s, Eq. (1)
   _specific_humidity_ratio[_qp] = _property_uo.specificHumidityRatio(_temperature[_qp]);
 
-  // _rho_{vs}, Eq. (3)
+  // P_{vs}, Eq. (2)
   _saturation_pressure_of_water_vapor_over_ice[_qp] = _property_uo.saturationPressureOfWaterVaporOverIce(_temperature[_qp]);
 }

src/userobjects/PropertyUserObject.C

@@ -59,7 +59,7 @@ PropertyUserObject::objectParams()
   params.addParam<Real>("atmospheric_pressure", 1.01325e5, "Atmospheric pressure, P_a [Pa]");
   params.addParam<Real>("gas_constant_dry_air", 286.9, "Gas constant for dry air, R_{da} [J/(Kg K)]");
   params.addParam<Real>("gas_constant_water_vapor", 461.5, "Gas constant for water vapor, R_v [J/(Kg K)]");
-  params.addParam<Real>("reference_temperature", 263.15 ,"Reference temperature, T_0 [K]");
+  params.addParam<Real>("reference_temperature", 263.15, "Reference temperature, T_0 [K]");
   params.addParam<Real>("interface_free_energy", 1.09e-1, "Interface free energy, \gamma [J/m^2]");
   params.addParam<Real>("mean_molecular_spacing", 3.19e-10, "Mean inter-molecular spacing in ice, a [m]");
   params.addParam<Real>("boltzmann", 1.3806488e-23, "Boltzmann constant, k [J/k]");
@@ -138,4 +138,3 @@ PropertyUserObject::equilibriumConcentration(const Real & T) const
   Real rho_vs_T_0 = equilibriumWaterVaporConcentrationAtSaturation(_T_0);
   return (rho_vs_T - rho_vs_T_0) / _rho_i;
 }
-

tests/materials/equilibrium_concentration/equilibrium_concentration.i

@@ -0,0 +1,60 @@
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+[]
+
+[Variables]
+  [./T]
+    initial_condition = 268.15
+  [../]
+[]
+
+[PikaMaterials]
+    conductivity_air = 0.1
+    temperature = T
+    phi = -1
+    output_properties = 'equilibrium_concentration density_air density_ice'
+    outputs = all
+  [../]
+[]
+
+[Postprocessors]
+  [./T]
+    type = PointValue
+    variable = T
+    point = '0.5 0.5 0'
+  [../]
+  [./u_eq]
+    type = PointValue
+    variable = 'equilibrium_concentration'
+    point = '0.5 0.5 0'
+  [../]
+[]
+
+[Problem]
+  type = FEProblem
+  solve = false
+  kernel_coverage_check = false
+[]
+
+[Executioner]
+  # Preconditioned JFNK (default)
+  type = Transient
+  num_steps = 1
+  solve_type = PJFNK
+  petsc_options_iname = '-pc_type -pc_hypre_type'
+  petsc_options_value = 'hypre boomeramg'
+[]
+
+[Outputs]
+  exodus = true
+  [./console]
+    type = Console
+    perf_log = true
+    nonlinear_residuals = true
+    linear_residuals = true
+  [../]
+  [./data]
+    type = CSV
+  [../]
+[]

tests/materials/equilibrium_concentration/gold/equilibrium_concentration_data.csv

@@ -0,0 +1,2 @@
+time,T,u_eq
+1,268.15,1.27847699784139e-6

tests/materials/equilibrium_concentration/tests

@@ -0,0 +1,8 @@
+[Tests]
+  [./equilibrium_concentration]
+    # Gold for this test comes python/KaempferPlapp2009.py
+    type = 'CSVDiff'
+    input = 'equilibrium_concentration.i'
+    csvdiff = 'equilibrium_concentration_data.csv'
+  [../]
+[]