From 000727f1e89b629ca81296e6dd51d596a2937ea4 Mon Sep 17 00:00:00 2001 From: Andrew E Slaughter Date: Wed, 18 Jun 2014 21:13:56 -0600 Subject: [PATCH] Added test for u_eq --- python/KaempferPlapp2009.py | 16 ++++- src/materials/PhaseFieldProperties.C | 3 +- src/userobjects/PropertyUserObject.C | 2 +- .../equilibrium_concentration.i | 60 +++++++++++++++++++ .../gold/equilibrium_concentration_data.csv | 2 + .../materials/equilibrium_concentration/tests | 8 +++ 6 files changed, 87 insertions(+), 4 deletions(-) create mode 100644 tests/materials/equilibrium_concentration/equilibrium_concentration.i create mode 100644 tests/materials/equilibrium_concentration/gold/equilibrium_concentration_data.csv create mode 100644 tests/materials/equilibrium_concentration/tests diff --git a/python/KaempferPlapp2009.py b/python/KaempferPlapp2009.py index 03d7b509..33e80691 100755 --- a/python/KaempferPlapp2009.py +++ b/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}) diff --git a/src/materials/PhaseFieldProperties.C b/src/materials/PhaseFieldProperties.C index 5682a28f..a3806658 100644 --- a/src/materials/PhaseFieldProperties.C +++ b/src/materials/PhaseFieldProperties.C @@ -57,7 +57,6 @@ PhaseFieldProperties::computeQpProperties() MaterialProperty & rho_i = getMaterialProperty("density_ice"); - /// @todo{This needs to be computed} _interface_velocity[_qp] = 1e-9; // [m/s] @@ -82,6 +81,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]); } diff --git a/src/userobjects/PropertyUserObject.C b/src/userobjects/PropertyUserObject.C index 12451db8..cd7c50f1 100644 --- a/src/userobjects/PropertyUserObject.C +++ b/src/userobjects/PropertyUserObject.C @@ -58,7 +58,7 @@ PropertyUserObject::objectParams() params.addParam("atmospheric_pressure", 1.01325e5, "Atmospheric pressure, P_a [Pa]"); params.addParam("gas_constant_dry_air", 286.9, "Gas constant for dry air, R_{da} [J/(Kg K)]"); params.addParam("gas_constant_water_vapor", 461.5, "Gas constant for water vapor, R_v [J/(Kg K)]"); - params.addParam("reference_temperature", 258.2 ,"Reference temperature, T_0 [K]"); + params.addParam("reference_temperature", 263.15,"Reference temperature, T_0 [K]"); params.addParam("interface_free_energy", 1.09e-1, "Interface free energy, \gamma [J/m^2]"); params.addParam("mean_molecular_spacing", 3.19e-10, "Mean inter-molecular spacing in ice, a [m]"); params.addParam("boltzmann", 1.3806488e-23, "Boltzmann constant, k [J/k]"); diff --git a/tests/materials/equilibrium_concentration/equilibrium_concentration.i b/tests/materials/equilibrium_concentration/equilibrium_concentration.i new file mode 100644 index 00000000..48d48ae0 --- /dev/null +++ b/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 + [../] +[] diff --git a/tests/materials/equilibrium_concentration/gold/equilibrium_concentration_data.csv b/tests/materials/equilibrium_concentration/gold/equilibrium_concentration_data.csv new file mode 100644 index 00000000..4a7eb61d --- /dev/null +++ b/tests/materials/equilibrium_concentration/gold/equilibrium_concentration_data.csv @@ -0,0 +1,2 @@ +time,T,u_eq +1,268.15,1.27847699784139e-6 diff --git a/tests/materials/equilibrium_concentration/tests b/tests/materials/equilibrium_concentration/tests new file mode 100644 index 00000000..c1eb4481 --- /dev/null +++ b/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' + [../] +[]