diff --git a/include/components/CoaxialPipe1Phase.h b/include/components/CoaxialPipe1Phase.h index bd1b644..863fa4f 100644 --- a/include/components/CoaxialPipe1Phase.h +++ b/include/components/CoaxialPipe1Phase.h @@ -23,8 +23,19 @@ class CoaxialPipe1Phase : public Coaxial1PhaseBase { // Add solid shell around annulus void AddSolidShell(const InputParameters ¶ms); + void addMooseObjects() override; + // Add solid-fluid connection based on component names void AddHeatTransferConnection(const std::string &flow_channel, const std::string &hs, const std::string &hs_side, const Real radius); -}; \ No newline at end of file + + // Add ambient convection to shell surface + void AddAmbientConvection(); + + // Create constant function based on scalar value + FunctionName CreateFunctionFromValue(const std::string &suffix, + const Real value); + + const Real _T_ambient, _p_ambient; +}; diff --git a/src/components/CoaxialPipe1Phase.C b/src/components/CoaxialPipe1Phase.C index b694d80..8286b43 100644 --- a/src/components/CoaxialPipe1Phase.C +++ b/src/components/CoaxialPipe1Phase.C @@ -1,12 +1,17 @@ #include "CoaxialPipe1Phase.h" #include "Component1D.h" +#include "Conversion.h" #include "FEProblemBase.h" #include "Factory.h" +#include "FlowModel.h" #include "InputParameters.h" +#include "MooseEnum.h" #include "MooseTypes.h" #include "Registry.h" +#include "SubProblem.h" #include "THMProblem.h" #include +#include registerMooseObject("ProteusApp", CoaxialPipe1Phase); @@ -101,6 +106,17 @@ InputParameters CoaxialPipe1Phase::validParams() { params.addParam( "outer_shell_Hw", "Manually specified HTC for annular pipe to shell."); + // Parameters for ambient convection + params.addParam("use_ambient_convection", false, + "Whether to apply ambient convection to the external " + "surface of the shell"); + params.addParam("T_ambient", 298, "Ambient temperature [K]."); + params.addParam("p_ambient", 101325, "Ambient pressure [Pa]."); + + MooseEnum ambient_properties("air", "air"); + params.addParam("ambient_properties", ambient_properties, + "Ambient fluid properties"); + // Add global parameter options params.addParam( "fp", "Global fluid properties. Overriden by inner_fp and outer_fp."); @@ -122,7 +138,8 @@ InputParameters CoaxialPipe1Phase::validParams() { } CoaxialPipe1Phase::CoaxialPipe1Phase(const InputParameters ¶ms) - : Coaxial1PhaseBase(params) { + : Coaxial1PhaseBase(params), _T_ambient(getParam("T_ambient")), + _p_ambient(getParam("p_ambient")) { // Add components AddInnerPipe(params); AddOuterAnnulus(params); @@ -140,6 +157,10 @@ CoaxialPipe1Phase::CoaxialPipe1Phase(const InputParameters ¶ms) AddHeatTransferConnection("outer", "tube", "OUTER", outer_radius); AddHeatTransferConnection("outer", "shell", "INNER", params.get("shell_inner_radius")); + + if (getParam("use_ambient_convection")) { + AddAmbientConvection(); + } } void CoaxialPipe1Phase::AddInnerPipe(const InputParameters ¶ms) { @@ -323,3 +344,144 @@ void CoaxialPipe1Phase::AddHeatTransferConnection( getTHMProblem().addComponent( class_name, name() + "_" + flow_channel + "_" + hs, ht_params); } + +void CoaxialPipe1Phase::AddAmbientConvection() { + + // Ambient convection + { + const std::string class_name = "HSBoundaryAmbientConvection"; + auto params = _factory.getValidParams(class_name); + params.set("_thm_problem") = &getTHMProblem(); + + params.set("T_ambient") = + CreateFunctionFromValue("T_ambient", _T_ambient); + params.set>("boundary") = {name() + + "/shell:outer"}; + params.set("hs") = {name() + "/shell"}; + params.set("htc_ambient") = "Hw"; + + getTHMProblem().addComponent(class_name, name() + "/conv_ambient", params); + } +} + +void CoaxialPipe1Phase::addMooseObjects() { + if (!getParam("use_ambient_convection")) + return; + + auto gravity = getParam("gravity_vector"); + if (gravity.norm() < 1e-8) + mooseError("The vertical vector must have magnitude greater than 0"); + + auto v = getParam("orientation"); + Real l{0.}; + Real dot_prod{fabs(v * gravity / (v.norm() * gravity.norm()))}; + std::string expression; + + if (dot_prod < 1e-8) { // horizontal pipe + auto widths = getParam>("shell_widths"); + l = 2. * std::accumulate(widths.begin(), widths.end(), + getParam("shell_inner_radius")); + expression = "pow(0.6 + (0.387*pow(Ra,1./6.))/pow(1 + " + "pow(0.559/Pr,9./16.),8./27.), 2)"; + + } else if (fabs(dot_prod - 1) < 1e-8) { // vertical pipe + auto lengths = getParam>("length"); + l = std::accumulate(lengths.begin(), lengths.end(), 0.); + expression = "pow(0.825 + (0.387*pow(Ra,1./6.))/pow(1 + " + "pow(0.492/Pr,9./16.),8/27), 2)"; + } else { + mooseError("Ambient convection only usable for vertical " + "and horizontal pipes."); + return; + } + + Real mu, k, cp, rho, R, gamma, beta; + if (getParam("ambient_properties") == "air") { + mu = 1.823e-05; + k = 0.02568; + R = 8.31446261815324 / 0.0289647; + gamma = 1.4; + cp = gamma * R / (gamma - 1); + rho = _p_ambient / (R * _T_ambient); + beta = 1 / _T_ambient; + } + + mooseInfo("Ambient properties\n", "\tk: ", k, "\n\tmu: ", mu, "\n\tcp: ", cp, + "\n\tbeta: ", beta, "\n\trho: ", rho, "\n"); + + // Create Rayleigh number property + { + const std::string class_name = "ADParsedFunctorMaterial"; + auto params = _factory.getValidParams(class_name); + params.set("_fe_problem_base") = &getTHMProblem(); + params.set("property_name") = "Ra"; + params.set("expression") = + "rho*beta*abs(T_solid-T_a)*L*L*L*g/(mu*k/(rho*cp))"; + params.set>("functor_symbols") = { + "rho", "beta", "mu", "k", "cp", "T_solid", "T_a", "L", "g"}; + + std::vector functor_names{ + Moose::stringifyExact(rho), + Moose::stringifyExact(beta), + Moose::stringifyExact(mu), + Moose::stringifyExact(k), + Moose::stringifyExact(cp), + "T_solid", + Moose::stringifyExact(_T_ambient), + Moose::stringifyExact(l), + Moose::stringifyExact(gravity.norm())}; + params.set>("functor_names") = functor_names; + params.set>("block") = { + name() + + "/shell:" + getParam>("shell_names").back()}; + + getTHMProblem().addMaterial(class_name, name() + "/Ra_conv", params); + } + + // Nusselt number + { + const std::string class_name = "ADParsedFunctorMaterial"; + auto params = _factory.getValidParams(class_name); + params.set("_fe_problem_base") = &getTHMProblem(); + params.set("property_name") = "Nu"; + + params.set("expression") = expression; + + params.set>("functor_symbols") = {"Pr", "Ra"}; + params.set>("functor_names") = { + Moose::stringifyExact(mu * cp / k), "Ra"}; + + params.set>("block") = { + name() + + "/shell:" + getParam>("shell_names").back()}; + getTHMProblem().addMaterial(class_name, name() + "/Nu_conv", params); + } + + // HTC + { + const std::string class_name = "ADParsedFunctorMaterial"; + auto params = _factory.getValidParams(class_name); + params.set("_fe_problem_base") = &getTHMProblem(); + params.set("property_name") = "Hw"; + params.set>("block") = { + name() + + "/shell:" + getParam>("shell_names").back()}; + + params.set("expression") = "Nu*k/L"; + params.set>("functor_symbols") = {"k", "L", "Nu"}; + params.set>("functor_names") = { + Moose::stringifyExact(k), Moose::stringifyExact(l), "Nu"}; + getTHMProblem().addMaterial(class_name, name() + "/Hw_conv", params); + } +} + +FunctionName +CoaxialPipe1Phase::CreateFunctionFromValue(const std::string &suffix, + const Real value) { + auto func_params = _factory.getValidParams("ConstantFunction"); + func_params.set("value") = value; + + auto func_name = name() + "_" + suffix; + getTHMProblem().addFunction("ConstantFunction", func_name, func_params); + return func_name; +} diff --git a/test/tests/components/coaxial_pipe/ambient_convection.i b/test/tests/components/coaxial_pipe/ambient_convection.i new file mode 100644 index 0000000..50a5668 --- /dev/null +++ b/test/tests/components/coaxial_pipe/ambient_convection.i @@ -0,0 +1,119 @@ +# Ambient convection correlation test + +T_solid = 350 +T_ambient = 300 + +[GlobalParams] + initial_p = 1e5 + initial_T = ${T_solid} + initial_vel = 0.1 + closures = thm_closures + fp = fluid +[] + +[FluidProperties] + [fluid] + type = SimpleFluidProperties + cv = 4000 + [] +[] + +[SolidProperties] + [solid] + type = ThermalFunctionSolidProperties + cp = 500 + k = 10 + rho = 1000 + [] +[] + +[Closures] + [thm_closures] + type = Closures1PhaseTHM + [] +[] + +[Components] + [inlet_inner] + type = InletMassFlowRateTemperature1Phase + T = ${T_solid} + m_dot = 0.1 + input = coaxial/inner:in + [] + [inlet_outer] + type = InletMassFlowRateTemperature1Phase + T = ${T_solid} + m_dot = 0.1 + input = coaxial/outer:in + [] + [coaxial] + type = CoaxialPipe1Phase + position = '0 0 0' + orientation = '1 0 0' + length = '0.4 0.6' + n_elems = '1 1' + axial_region_names = 'section_1 section_2' + + tube_inner_radius = 0.025 + tube_names = tube + tube_widths = 0.025 + tube_materials = solid + tube_n_elems = 1 + tube_T_ref = ${T_solid} + + shell_inner_radius = 0.075 + shell_names = shell + shell_widths = 0.025 + shell_materials = solid + shell_n_elems = 1 + shell_T_ref = ${T_solid} + + use_ambient_convection = true + T_ambient = ${T_ambient} + [] + [outlet_inner] + type = Outlet1Phase + input = coaxial/inner:out + p = 1e5 + [] + [outlet_outer] + type = Outlet1Phase + input = coaxial/outer:out + p = 1e5 + [] +[] + +[Postprocessors] + [Ra] + type = ADElementExtremeFunctorValue + functor = Ra + block = coaxial/shell:shell + execute_on = INITIAL + [] + [Nu] + type = ADElementExtremeFunctorValue + functor = Nu + block = coaxial/shell:shell + execute_on = INITIAL + [] + [Hw] + type = ADElementExtremeFunctorValue + functor = Hw + block = coaxial/shell:shell + execute_on = INITIAL + [] +[] + +[Problem] + solve = false +[] + +[Executioner] + type = Steady +[] + +[Outputs] + csv = true + execute_on = INITIAL + show = 'Ra Nu Hw' +[] diff --git a/test/tests/components/coaxial_pipe/test.py b/test/tests/components/coaxial_pipe/test.py index c21470a..4ca3675 100644 --- a/test/tests/components/coaxial_pipe/test.py +++ b/test/tests/components/coaxial_pipe/test.py @@ -3,6 +3,34 @@ import unittest import numpy as np + +def ambient_convection_values(T_solid, T_ambient, length, gravity, vertical): + """Independently evaluate the ambient convection correlations.""" + mu = 1.823e-5 + k = 0.02568 + gas_constant = 8.31446261815324 / 0.0289647 + gamma = 1.4 + cp = gamma * gas_constant / (gamma - 1) + rho = 101325 / (gas_constant * T_ambient) + beta = 1 / T_ambient + + prandtl = mu * cp / k + thermal_diffusivity = k / (rho * cp) + rayleigh = (rho * beta * abs(T_solid - T_ambient) * length**3 * gravity + / (mu * thermal_diffusivity)) + + base, prandtl_coefficient = (0.825, 0.492) if vertical else (0.6, 0.559) + nusselt = (base + 0.387 * rayleigh ** (1 / 6) + / (1 + (prandtl_coefficient / prandtl) ** (9 / 16)) ** (8 / 27)) ** 2 + + return np.array([nusselt * k / length, nusselt, rayleigh]) + + +def read_ambient_convection_output(file_name): + """Read Hw, Nu, and Ra from an ambient convection CSV output.""" + data = np.genfromtxt(file_name, delimiter=",", names=True) + return np.array([data["Hw"], data["Nu"], data["Ra"]]) + class TestCoaxialPipe(unittest.TestCase): """Test class for the coaxial pipe component.""" def test_energy_balance(self): @@ -74,3 +102,43 @@ def test_energy_balance_outer(self): rel_diff = abs(total_energy - q)/q assert rel_diff < 0.00048, f"Rel. energy difference greater than 0.00046: {rel_diff}" + + def assert_ambient_convection_values( + self, file_name, T_solid, T_ambient, length, gravity, vertical + ): + """Compare computed correlation values against an independent evaluation.""" + actual = read_ambient_convection_output(file_name) + expected = ambient_convection_values( + T_solid, T_ambient, length, gravity, vertical + ) + np.testing.assert_allclose(actual, expected, rtol=1e-12) + + def test_ambient_convection_horizontal(self): + """Checks the diameter-based horizontal cylinder correlation.""" + self.assert_ambient_convection_values( + "ambient_convection_horizontal.csv", 350, 300, 0.2, 9.81, False + ) + + def test_ambient_convection_vertical(self): + """Checks the total-length-based vertical correlation.""" + self.assert_ambient_convection_values( + "ambient_convection_vertical.csv", 350, 300, 1.0, 9.81, True + ) + + def test_ambient_convection_vertical_reversed(self): + """Checks that reversing the vertical direction does not change the result.""" + self.assert_ambient_convection_values( + "ambient_convection_vertical_reversed.csv", 350, 300, 1.0, 9.81, True + ) + + def test_ambient_convection_cold_surface(self): + """Checks a pipe colder than its ambient environment.""" + self.assert_ambient_convection_values( + "ambient_convection_cold_surface.csv", 250, 300, 0.2, 9.81, False + ) + + def test_ambient_convection_gravity_magnitude(self): + """Checks that the configured gravity magnitude is used.""" + self.assert_ambient_convection_values( + "ambient_convection_gravity_magnitude.csv", 350, 300, 0.2, 4.905, False + ) diff --git a/test/tests/components/coaxial_pipe/tests b/test/tests/components/coaxial_pipe/tests index f6b8225..35f945c 100644 --- a/test/tests/components/coaxial_pipe/tests +++ b/test/tests/components/coaxial_pipe/tests @@ -31,7 +31,7 @@ [test] type= PythonUnitTest input = test.py - test_case = TestCoaxialPipe.test_energy_balance_outer + test_case = TestCoaxialPipe.test_energy_balance_inner heavy=true prereq = energy_balance_inner/run [] @@ -50,4 +50,90 @@ [] requirement = "Checks that outer pipe conserved energy in isolation" [] + [ambient_convection] + [horizontal] + [run] + type = RunApp + input = ambient_convection.i + cli_args = "Outputs/file_base=ambient_convection_horizontal" + [] + [test] + type = PythonUnitTest + input = test.py + test_case = TestCoaxialPipe.test_ambient_convection_horizontal + prereq = ambient_convection/horizontal/run + [] + requirement = "The system shall compute ambient natural convection for a horizontal coaxial pipe using its outer diameter." + [] + [vertical] + [run] + type = RunApp + input = ambient_convection.i + cli_args = "Components/coaxial/orientation='0 0 1' Outputs/file_base=ambient_convection_vertical" + [] + [test] + type = PythonUnitTest + input = test.py + test_case = TestCoaxialPipe.test_ambient_convection_vertical + prereq = ambient_convection/vertical/run + [] + requirement = "The system shall compute ambient natural convection for a vertical coaxial pipe using its total axial length." + [] + [vertical_reversed] + [run] + type = RunApp + input = ambient_convection.i + cli_args = "Components/coaxial/orientation='0 0 -1' Outputs/file_base=ambient_convection_vertical_reversed" + [] + [test] + type = PythonUnitTest + input = test.py + test_case = TestCoaxialPipe.test_ambient_convection_vertical_reversed + prereq = ambient_convection/vertical_reversed/run + [] + requirement = "The system shall compute identical ambient convection for either vertical pipe direction." + [] + [cold_surface] + [run] + type = RunApp + input = ambient_convection.i + cli_args = "T_solid=250 Outputs/file_base=ambient_convection_cold_surface" + [] + [test] + type = PythonUnitTest + input = test.py + test_case = TestCoaxialPipe.test_ambient_convection_cold_surface + prereq = ambient_convection/cold_surface/run + [] + requirement = "The system shall compute finite ambient convection properties when the pipe is colder than its environment." + [] + [gravity_magnitude] + [run] + type = RunApp + input = ambient_convection.i + cli_args = "Components/coaxial/gravity_vector='0 0 -4.905' Outputs/file_base=ambient_convection_gravity_magnitude" + [] + [test] + type = PythonUnitTest + input = test.py + test_case = TestCoaxialPipe.test_ambient_convection_gravity_magnitude + prereq = ambient_convection/gravity_magnitude/run + [] + requirement = "The system shall use the configured gravity magnitude in the ambient convection correlation." + [] + [oblique] + type = RunException + input = ambient_convection.i + cli_args = "Components/coaxial/orientation='1 0 1'" + expect_err = "Ambient convection only usable for vertical and horizontal pipes." + requirement = "The system shall reject ambient convection for an oblique coaxial pipe." + [] + [zero_gravity] + type = RunException + input = ambient_convection.i + cli_args = "Components/coaxial/gravity_vector='0 0 0'" + expect_err = "magnitude greater than 0" + requirement = "The system shall reject ambient convection when gravity has zero magnitude." + [] + [] []