Radial Heat Conduction Verification

Test Description

This verification case checks the radial heat conduction between subchannels in the triangular-lattice solver, TriSubChannel1PhaseProblem. In the enthalpy equation of SubChannel Theory each gap $var element = document.getElementById("moose-equation-2b388ec0-1ae3-4a3f-8a3d-2765ca884e80");katex.render("ij", element, {displayMode:false,throwOnError:false});$ between subchannels $var element = document.getElementById("moose-equation-b64d9dc4-8617-40c4-8d13-57a4e7ec8ea5");katex.render("i", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-db8464a2-93b1-423e-b9cb-d4fdc1d75198");katex.render("j", element, {displayMode:false,throwOnError:false});$ conducts

$var element = document.getElementById("moose-equation-fce63851-3f45-4382-ae17-2ab2e830a27e");katex.render("Q_{ij,k} = k \\frac{S_{ij,k} \\eta_{ij}}{L_{ij}} (T_{i,k} - T_{j,k}),", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-ddd47867-8c4b-45c3-b66e-71e821e99829");katex.render("k", element, {displayMode:false,throwOnError:false});$ is the fluid thermal conductivity and $var element = document.getElementById("moose-equation-50c4acaf-e573-4ce5-a3ca-a9d4c7ecb42a");katex.render("S_{ij,k} = w_{ij} \\Delta z_k", element, {displayMode:false,throwOnError:false});$ is the gap area. The gap width $var element = document.getElementById("moose-equation-d6892369-9911-4fef-b944-871f99f3a948");katex.render("w_{ij}", element, {displayMode:false,throwOnError:false});$ is $var element = document.getElementById("moose-equation-52bd60f2-6b70-4d18-adef-5a948ba2f823");katex.render("p - D", element, {displayMode:false,throwOnError:false});$ between two pins and the duct-to-pin gap where the gap faces the duct, with pitch $var element = document.getElementById("moose-equation-4a8b87a6-9948-4d68-8885-b0e166cc3cc7");katex.render("p", element, {displayMode:false,throwOnError:false});$ and pin diameter $var element = document.getElementById("moose-equation-ad7a6e98-0299-445a-8202-c0272ace4d89");katex.render("D", element, {displayMode:false,throwOnError:false});$ . The centroid distance $var element = document.getElementById("moose-equation-9d303b1f-63bc-4d9d-8ea6-854c883b21dc");katex.render("L_{ij}", element, {displayMode:false,throwOnError:false});$ is $var element = document.getElementById("moose-equation-cdb8cfee-2208-4613-9a92-051c89fe8d93");katex.render("p / \\sqrt{3}", element, {displayMode:false,throwOnError:false});$ for gaps that involve an interior subchannel and $var element = document.getElementById("moose-equation-8bc28f2e-db24-4f19-ac25-4cc7eeaa1983");katex.render("p", element, {displayMode:false,throwOnError:false});$ otherwise, and the shape factor is

$var element = document.getElementById("moose-equation-d8a97089-1939-4869-b5ab-f515e4ee05ae");katex.render("\\eta_{ij} = 0.66 \\frac{p}{D} \\left( \\frac{w_{ij}}{D} \\right)^{-0.3}.", element, {displayMode:true,throwOnError:false});$

The problem isolates the conduction term. It consists of a 7-pin assembly with $var element = document.getElementById("moose-equation-1b39c686-c38f-4a39-8a6b-d27517155b9a");katex.render("p = 12", element, {displayMode:false,throwOnError:false});$ mm, $var element = document.getElementById("moose-equation-4aa9f651-d9df-48f6-8561-0a301b039d34");katex.render("D = 10", element, {displayMode:false,throwOnError:false});$ mm, a duct flat-to-flat distance of 35 mm, and a heated length of 0.5 m divided into 20 axial cells. The pins are unheated and bare. Turbulent mixing is disabled ( $var element = document.getElementById("moose-equation-58c6d44d-767d-49e9-a4f5-93806d96c52d");katex.render("\\beta = 0", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-b37af5c3-332e-4f09-ad9c-e94da74ece17");katex.render("C_T = 0", element, {displayMode:false,throwOnError:false});$ ), and a lateral form loss coefficient $var element = document.getElementById("moose-equation-12d0a4b4-cdcd-4ea7-b428-7d8fc5c08d88");katex.render("K_{ij} = 10^8", element, {displayMode:false,throwOnError:false});$ suppresses diversion crossflow. The fluid properties are constant and close to those of sodium at 660 K: $var element = document.getElementById("moose-equation-737687ab-fbae-4362-974e-733ef20195ae");katex.render("\\rho = 850", element, {displayMode:false,throwOnError:false});$ kg/m $var element = document.getElementById("moose-equation-b36fd73f-aca9-499e-9405-6033006e225e");katex.render("^3", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-b33476d0-c94d-44c4-859f-25ea0792c654");katex.render("c_p = 1270", element, {displayMode:false,throwOnError:false});$ J/kg-K, $var element = document.getElementById("moose-equation-af9ea54c-7c0d-475c-afb8-3d130b4b3199");katex.render("k = 70", element, {displayMode:false,throwOnError:false});$ W/m-K, and $var element = document.getElementById("moose-equation-e134964b-64a2-4ff8-9b46-fda60ed0efdf");katex.render("\\mu = 2.5 \\times 10^{-4}", element, {displayMode:false,throwOnError:false});$ Pa-s, so $var element = document.getElementById("moose-equation-a2090a37-5edb-4e47-b3c4-c3a60d0795ba");katex.render("h = c_p T", element, {displayMode:false,throwOnError:false});$ . The duct heats the subchannels with a heat flux of $var element = document.getElementById("moose-equation-9185f1b0-bc5f-44be-8425-0f8e1c2796e9");katex.render("2 \\times 10^4", element, {displayMode:false,throwOnError:false});$ W/m $var element = document.getElementById("moose-equation-0ca47d10-4c5f-4596-a954-93122ec65a0f");katex.render("^2", element, {displayMode:false,throwOnError:false});$ where $var element = document.getElementById("moose-equation-2193d22e-e9e8-4fad-ae16-153f0ba42f83");katex.render("x > 0", element, {displayMode:false,throwOnError:false});$ and cools them with the same flux where $var element = document.getElementById("moose-equation-422684b7-61f5-41ac-9fad-3007b549143c");katex.render("x < 0", element, {displayMode:false,throwOnError:false});$ , so the assembly receives no net heat and heat crosses it only by conduction through the gaps. The inlet mass flux is 5 kg/m $var element = document.getElementById("moose-equation-a129be95-6cf0-43d1-84ae-d0c2522631f6");katex.render("^2", element, {displayMode:false,throwOnError:false});$ -s, low enough that the outlet temperatures span about 50 K.

The reference solution solves the same conductance network independently of SCM. It rebuilds the 24 gaps of the assembly from the subchannel positions, assembles the conductances above, adds the duct heat of each edge and corner subchannel as the heat flux times its peripheral duct width, and marches the enthalpy cell by cell, implicitly in the conduction term. It takes only the subchannel positions, flow areas, and face mass flow rates from SCM, and not the SCM matrices. The flow split develops over the first cells, where the axial mass flow rate redistributes between subchannels; the reference solution uses the SCM face mass flow rates but does not model the enthalpy carried by that redistribution.

Results

Figure 1 compares the outlet subchannel temperatures. Heat flows from the heated corner subchannel at $var element = document.getElementById("moose-equation-693def00-969b-498e-9b24-be326da1e6fa");katex.render("x = 18", element, {displayMode:false,throwOnError:false});$ mm, 24.8 K above the inlet temperature, to the cooled corner subchannel at $var element = document.getElementById("moose-equation-519635d0-1928-4010-a731-69519ed92d9b");katex.render("x = -18", element, {displayMode:false,throwOnError:false});$ mm, 24.8 K below it, and the subchannels on $var element = document.getElementById("moose-equation-602427a1-f6da-4079-a6f4-6f45c0307fda");katex.render("x = 0", element, {displayMode:false,throwOnError:false});$ stay at the inlet temperature. SCM and the network solution agree within $var element = document.getElementById("moose-equation-592861ae-d15d-4cbe-aa02-2e798b3443a6");katex.render("3.4 \\times 10^{-3}", element, {displayMode:false,throwOnError:false});$ K, about $var element = document.getElementById("moose-equation-52da0ba3-ec6c-4e3d-b702-a3307d9e14ec");katex.render("10^{-4}", element, {displayMode:false,throwOnError:false});$ of the temperature spread. The remaining difference is largest in the interior subchannels and comes from the flow redistribution over the development length.

Left: outlet subchannel temperature above the inlet temperature from SCM and from the independent network solution, versus the subchannel centroid position $x$. Right: difference between SCM and the network solution.

Figure 1: Left: outlet subchannel temperature above the inlet temperature from SCM and from the independent network solution, versus the subchannel centroid position $var element = document.getElementById("moose-equation-f396ebd9-0850-49a0-a4fd-b54d09b81609");katex.render("x", element, {displayMode:false,throwOnError:false});$ . Right: difference between SCM and the network solution.

Input file

The verification problem uses the following input file:

# Radial heat conduction verification in a 7-pin triangular assembly. The pins are unheated and
# bare, turbulent mixing is disabled (beta = 0, CT = 0), and a large lateral form loss
# coefficient Kij suppresses diversion crossflow. The duct heats the subchannels on one side of
# the assembly and cools those on the other side with the same total power, so heat crosses the
# assembly only by radial conduction through the gaps. The axial flow is low so that the
# temperature differences between subchannels are large. The outlet subchannel temperatures are
# compared with an independent solution of the same gap conductance network in
# test_radial_conduction.py.
T_in = 660 # K
mass_flux_in = 5.0 # kg/m2/s
P_out = 2.0e5 # Pa
duct_flux = 2.0e4 # W/m2

[GlobalParams<<<{"href": "../../../syntax/GlobalParams/index.html"}>>>]
  nrings = 2
  n_cells = 20
  flat_to_flat = 0.035
  heated_length = 0.5
  pitch = 0.012
[]

[TriSubChannelMesh<<<{"href": "../../../syntax/TriSubChannelMesh/index.html"}>>>]
  [subchannel]
    type = SCMTriAssemblyMeshGenerator<<<{"description": "Creates a mesh of 1D subchannels and 1D pins in a triangular lattice arrangement", "href": "../../../source/meshgenerators/SCMTriAssemblyMeshGenerator.html"}>>>
    pin_diameter<<<{"description": "Rod diameter [m]"}>>> = 0.01
    dwire<<<{"description": "Wire diameter [m]"}>>> = 0.0
    hwire<<<{"description": "Wire lead length [m]"}>>> = 0.0
    spacer_z<<<{"description": "Axial location of spacers/vanes/mixing vanes [m]"}>>> = '0'
    spacer_k<<<{"description": "K-loss coefficient of spacers/vanes/mixing vanes [-]"}>>> = '0'
    Kij<<<{"description": "Lateral form loss coefficient [-]"}>>> = 1e8
  []
  [duct]
    type = SCMTriDuctMeshGenerator<<<{"description": "Creates a mesh of 2D duct cells around a triangular lattice subassembly", "href": "../../../source/meshgenerators/SCMTriDuctMeshGenerator.html"}>>>
    input<<<{"description": "The corresponding subchannel mesh"}>>> = subchannel
  []
[]

# Constant properties close to sodium at 660 K, so that the network solution has constant
# conductances. With cp = cv, zero thermal expansion and porepressure_coefficient = 0,
# h = cp * T exactly
[FluidProperties<<<{"href": "../../../syntax/FluidProperties/index.html"}>>>]
  [sodium]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 850
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
    cp<<<{"description": "Constant specific heat capacity at constant pressure (J/kg/K)"}>>> = 1270
    cv<<<{"description": "Constant specific heat capacity at constant volume (J/kg/K)"}>>> = 1270
    thermal_conductivity<<<{"description": "Constant thermal conductivity (W/m/K)"}>>> = 70
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 2.5e-4
    porepressure_coefficient<<<{"description": "The enthalpy is internal_energy + P / density * porepressure_coefficient.  Physically this should be 1.0, but analytic solutions are simplified when it is zero"}>>> = 0
  []
[]

[SubChannel<<<{"href": "../../../syntax/SubChannel/index.html"}>>>]
  type = TriSubChannel1PhaseProblem
  fp = sodium
  P_out = ${P_out}
  compute_density = false
  compute_viscosity = false
  compute_power = true
  implicit = true
  segregated = true
  rtol = 1e-10
  atol = 1e-10
  friction_closure = 'Chen'
  mixing_closure = 'no_mixing'
  # The pin heat transfer closure only sets the pin surface temperature after the solve; with
  # q_prime = 0 the pins add no heat to the coolant. It is required whenever a pin mesh exists
  pin_HTC_closure = 'gnielinski'
  duct_HTC_closure = 'gnielinski'
[]

[SCMClosures<<<{"href": "../../../syntax/SCMClosures/index.html"}>>>]
  [Chen]
    type = SCMFrictionChenTodreas<<<{"description": "Class that computes the axial friction factor using the Chen-Todreas correlations.", "href": "../../../source/scmclosures/SCMFrictionChenTodreas.html"}>>>
  []
  [gnielinski]
    type = SCMHTCGnielinski<<<{"description": "Class that computes the convective heat transfer coefficient using the Gnielinski correlation.", "href": "../../../source/scmclosures/SCMHTCGnielinski.html"}>>>
  []
  [no_mixing]
    type = SCMMixingConstantBeta<<<{"description": "Class that models the turbulent mixing coefficient beta as a user defined constant.", "href": "../../../source/scmclosures/SCMMixingConstantBeta.html"}>>>
    beta<<<{"description": "Turbulent mixing parameter [-]."}>>> = 0.0
    CT<<<{"description": "Turbulent momentum, modeling parameter [-]."}>>> = 0.0
  []
[]

[Functions<<<{"href": "../../../syntax/Functions/index.html"}>>>]
  # Heat the duct where x > 0 and cool it where x < 0. The duct nodes on x = 0 get no flux, so
  # the heated and cooled sides receive the same power
  [duct_flux_fn]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = 'if(x > 1e-6, ${duct_flux}, if(x < -1e-6, -${duct_flux}, 0))'
  []
[]

[ICs<<<{"href": "../../../syntax/Cardinal/ICs/index.html"}>>>]
  [T_ic]
    type = ConstantIC<<<{"description": "Sets a constant field value.", "href": "../../../source/ics/ConstantIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = T
    value<<<{"description": "The value to be set in IC"}>>> = ${T_in}
  []
  [P_ic]
    type = ConstantIC<<<{"description": "Sets a constant field value.", "href": "../../../source/ics/ConstantIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = P
    value<<<{"description": "The value to be set in IC"}>>> = 0.0
  []
  [mu_ic]
    type = ViscosityIC<<<{"description": "Computes viscosity from specified pressure and temperature", "href": "../../../source/ics/ViscosityIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = mu
    p<<<{"description": "Pressure [Pa]"}>>> = ${P_out}
    T<<<{"description": "Temperature [K]"}>>> = T
    fp<<<{"description": "Fluid properties user object name"}>>> = sodium
  []
  [rho_ic]
    type = RhoFromPressureTemperatureIC<<<{"description": "Computes the density from pressure and temperature.", "href": "../../../source/ics/RhoFromPressureTemperatureIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = rho
    p<<<{"description": "The pressure [Pa]"}>>> = ${P_out}
    T<<<{"description": "The temperature [K]"}>>> = T
    fp<<<{"description": "The name of fluid properties user object."}>>> = sodium
  []
  [h_ic]
    type = SpecificEnthalpyFromPressureTemperatureIC<<<{"description": "Computes the specific enthalpy from pressure and temperature.", "href": "../../../source/ics/SpecificEnthalpyFromPressureTemperatureIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = h
    p<<<{"description": "The pressure [Pa]"}>>> = ${P_out}
    T<<<{"description": "The temperature [K]"}>>> = T
    fp<<<{"description": "The name of fluid properties user object."}>>> = sodium
  []
  [mdot_ic]
    type = ConstantIC<<<{"description": "Sets a constant field value.", "href": "../../../source/ics/ConstantIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = mdot
    value<<<{"description": "The value to be set in IC"}>>> = 0.0
  []
  [q_prime_ic]
    type = ConstantIC<<<{"description": "Sets a constant field value.", "href": "../../../source/ics/ConstantIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = q_prime
    value<<<{"description": "The value to be set in IC"}>>> = 0.0
  []
  [Dpin_ic]
    type = ConstantIC<<<{"description": "Sets a constant field value.", "href": "../../../source/ics/ConstantIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = Dpin
    value<<<{"description": "The value to be set in IC"}>>> = 0.01
  []
  [duct_heat_flux_ic]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = duct_heat_flux
    function<<<{"description": "The initial condition function."}>>> = duct_flux_fn
  []
[]

[AuxKernels<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
  [T_in_bc]
    type = ConstantAux<<<{"description": "Creates a constant field in the domain.", "href": "../../../source/auxkernels/ConstantAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = T
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = inlet
    value<<<{"description": "Some constant value that can be read from the input file"}>>> = ${T_in}
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_begin'
  []
  [mdot_in_bc]
    type = SCMMassFlowRateAux<<<{"description": "Computes mass flow rate from specified mass flux and subchannel cross-sectional area. Can read either PostprocessorValue or Real", "href": "../../../source/auxkernels/SCMMassFlowRateAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = mdot
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = inlet
    area<<<{"description": "Cross sectional area [m^2]"}>>> = S
    mass_flux<<<{"description": "The postprocessor or Real to use for the value of mass_flux"}>>> = ${mass_flux_in}
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_begin'
  []
[]

[VectorPostprocessors<<<{"href": "../../../syntax/VectorPostprocessors/index.html"}>>>]
  # Sampled at every axial level for the independent network solution in
  # test_radial_conduction.py
  [channels]
    type = NodalValueSampler<<<{"description": "Samples values of nodal variable(s).", "href": "../../../source/vectorpostprocessors/NodalValueSampler.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = 'T mdot S'
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = subchannel
    sort_by<<<{"description": "What to sort the samples by. Options include 'x', 'y', 'z', 'id', and the name of any of the sampled quantities (which each create a vector of the same name)."}>>> = id
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [duct]
    type = NodalValueSampler<<<{"description": "Samples values of nodal variable(s).", "href": "../../../source/vectorpostprocessors/NodalValueSampler.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = 'duct_heat_flux'
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = duct
    sort_by<<<{"description": "What to sort the samples by. Options include 'x', 'y', 'z', 'id', and the name of any of the sampled quantities (which each create a vector of the same name)."}>>> = id
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]

[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [Total_power_SCMTHPowerPostprocessor]
    type = SCMTHPowerPostprocessor<<<{"description": "Calculates the total power of the subchannel assembly $[W]$ via the thermal-hydraulic enthalpy flow rate balance between inlet and outlet.", "href": "../../../source/postprocessors/SCMTHPowerPostprocessor.html"}>>>
  []
  [Total_power_SCMDuctPowerPostprocessor]
    type = SCMDuctHeatRatePostprocessor<<<{"description": "Calculates the heat flow that goes into the coolant through the duct based on aux variable duct_heat_flux", "href": "../../../source/postprocessors/SCMDuctHeatRatePostprocessor.html"}>>>
  []
[]

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Steady
[]

[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
[]

The test test_radial_conduction.py runs it, solves the conductance network, and checks that the outlet temperatures agree within $var element = document.getElementById("moose-equation-35c347cc-0b20-4fce-adca-eefd9ab22e00");katex.render("4 \\times 10^{-3}", element, {displayMode:false,throwOnError:false});$ K:

"""
Radial heat conduction verification for the triangular-lattice subchannel solver.

radial_conduction_7pin.i heats the duct on one side of a 7-pin assembly and cools it on the other
with the same power. The pins are unheated, turbulent mixing is off, the fluid properties are
constant, and h = cp * T, so heat crosses the assembly only through the gap conductances.

This test solves the same conductance network independently of the solver and compares the
outlet subchannel temperatures. The network uses the documented radial conduction model of
TriSubChannel1PhaseProblem: each gap between channels i and j carries
k * w_ij * F_ij * dz / d_ij * (T_i - T_j), with F_ij = 0.66 (p / D) (w_ij / D)^-0.3. The gap width
w_ij is p - D for gaps between pins and the duct-to-pin gap for gaps facing the duct, and the
centroid distance d_ij is p / sqrt(3) for gaps involving an interior channel and p otherwise. The
duct adds q'' times the peripheral duct width of each edge or corner channel. The network is
marched in z with the solver's face mass flow rates, implicitly in the conduction term, which is
the discrete model the solver assembles. Only the channel positions, flow areas and mass flow
rates are taken from the solver output; the gap connectivity is rebuilt from the lattice.
"""

import csv
import math
import os
import subprocess
import unittest

import numpy as np

class TestRadialConductionVerification(unittest.TestCase):

    # Problem parameters (must match radial_conduction_7pin.i)
    T_in = 660.0
    cp = 1270.0
    k = 70.0
    pitch = 0.012
    pin_diameter = 0.01
    flat_to_flat = 0.035
    nrings = 2

    @classmethod
    def setUpClass(cls):
        cls.test_dir = os.path.dirname(os.path.abspath(__file__))
        cls.app = os.environ.get(
            "MOOSE_PYTHONUNITTEST_EXECUTABLE",
            os.path.join(cls.test_dir, "..", "..", "subchannel-opt"),
        )
        cls.input = os.path.join(cls.test_dir, "radial_conduction_7pin.i")
        cls.duct_to_pin_gap = 0.5 * (
            cls.flat_to_flat
            - (cls.nrings - 1) * cls.pitch * math.sqrt(3.0)
            - cls.pin_diameter
        )

    def _run(self, outbase, extra_args=None):
        """Run the solver and return the channel and duct samples at every axial level."""
        cmd = [
            self.app,
            "-i",
            self.input,
            f"Outputs/file_base={outbase}",
            "--no-color",
        ] + (extra_args or [])
        proc = subprocess.run(cmd, capture_output=True, text=True, cwd=self.test_dir)
        if proc.returncode != 0:
            self.fail(
                f"Solver failed ({outbase}):\n"
                f"--- stdout ---\n{proc.stdout[-3000:]}\n"
                f"--- stderr ---\n{proc.stderr[-1000:]}"
            )
        with open(os.path.join(self.test_dir, f"{outbase}_channels_0001.csv")) as fh:
            channels = list(csv.DictReader(fh))
        with open(os.path.join(self.test_dir, f"{outbase}_duct_0001.csv")) as fh:
            duct = list(csv.DictReader(fh))
        return channels, duct

    def _network_outlet_temperature(self, channels, duct):
        """Solve the gap conductance network; return (positions, T_solver, T_network)."""
        p, D, w_duct = self.pitch, self.pin_diameter, self.duct_to_pin_gap

        nodes = {}
        for r in channels:
            key = (round(float(r["x"]), 6), round(float(r["y"]), 6))
            nodes.setdefault(key, {})[round(float(r["z"]), 9)] = r
        positions = sorted(nodes)
        zs = sorted(nodes[positions[0]])
        n = len(positions)
        xy = np.array(positions)

        # Interior channels lie inside the outer pin ring; of the peripheral channels, the corner
        # channels have the smallest flow area
        area = np.array([float(nodes[c][zs[0]]["S"]) for c in positions])
        radius = np.hypot(xy[:, 0], xy[:, 1])
        r_outer_pins = (self.nrings - 1) * p
        kind = [
            (
                "interior"
                if r < r_outer_pins
                else ("corner" if math.isclose(s, area.min(), rel_tol=1e-9) else "edge")
            )
            for r, s in zip(radius, area)
        ]

        # Two channels share a gap when their centroids are the closest pair of their kinds.
        # Interior channels connect to interior and edge channels, edge channels to corner
        # channels; in a 7-pin assembly no two edge channels share a gap
        allowed = {("interior", "interior"), ("edge", "interior"), ("corner", "edge")}
        pairs = {}
        for i in range(n):
            for j in range(i + 1, n):
                kinds = tuple(sorted((kind[i], kind[j])))
                if kinds in allowed:
                    pairs.setdefault(kinds, []).append((math.dist(xy[i], xy[j]), i, j))
        gaps = []
        for kinds, candidates in pairs.items():
            d_min = min(c[0] for c in candidates)
            gaps += [(i, j, kinds) for d, i, j in candidates if d < 1.05 * d_min]

        # Conductance per unit length of each gap
        G = np.zeros((n, n))
        for i, j, kinds in gaps:
            if "interior" in kinds:
                width, dist = p - D, p / math.sqrt(3.0)
            else:
                width, dist = w_duct, p
            shape_factor = 0.66 * (p / D) * (width / D) ** -0.3
            g = self.k * width * shape_factor / dist
            G[i, i] += g
            G[j, j] += g
            G[i, j] -= g
            G[j, i] -= g

        # Duct heat per unit length: flux at the nearest duct node times the peripheral width
        duct_by_z = {}
        for r in duct:
            duct_by_z.setdefault(round(float(r["z"]), 9), []).append(r)
        corner_width = 2.0 / math.sqrt(3.0) * (D / 2.0 + w_duct)

        def duct_heat(i, z):
            if kind[i] == "interior":
                return 0.0
            x, y = positions[i]
            nearest = min(
                duct_by_z[z],
                key=lambda r: math.hypot(float(r["x"]) - x, float(r["y"]) - y),
            )
            width = p if kind[i] == "edge" else corner_width
            return float(nearest["duct_heat_flux"]) * width

        def mdot(z):
            return np.array([float(nodes[c][z]["mdot"]) for c in positions])

        h = np.full(n, self.cp * self.T_in)
        for iz in range(1, len(zs)):
            dz = zs[iz] - zs[iz - 1]
            q = np.array(
                [
                    0.5 * (duct_heat(i, zs[iz - 1]) + duct_heat(i, zs[iz])) * dz
                    for i in range(n)
                ]
            )
            h = np.linalg.solve(
                np.diag(mdot(zs[iz])) + G * dz / self.cp, mdot(zs[iz - 1]) * h + q
            )

        T_solver = np.array([float(nodes[c][zs[-1]]["T"]) for c in positions])
        return positions, T_solver, h / self.cp

    def test_matches_conductance_network(self):
        """Outlet temperatures agree with the independent network solution."""
        positions, T_solver, T_network = self._network_outlet_temperature(
            *self._run("radial_conduction_7pin_out")
        )
        spread = np.abs(T_solver - self.T_in).max()
        error = np.abs(T_solver - T_network).max()

        print(f"\nRadial conduction verification, 7-pin assembly")
        print(f'{"x [m]":>9} {"y [m]":>9} {"T solver [K]":>13} {"T network [K]":>14}')
        for (x, y), a, b in sorted(zip(positions, T_solver, T_network)):
            print(f"{x:+9.5f} {y:+9.5f} {a:13.5f} {b:14.5f}")
        print(f"Largest |T - T_in|: {spread:.3f} K, largest difference: {error:.2e} K")

        # Write the comparison for the verification page
        with open(
            os.path.join(self.test_dir, "radial_conduction_7pin_comparison.csv"), "w"
        ) as fh:
            fh.write("x,y,T_solver,T_network\n")
            for (x, y), a, b in sorted(zip(positions, T_solver, T_network)):
                fh.write(f"{x:.6f},{y:.6f},{a:.6f},{b:.6f}\n")

        # The heated and cooled sides span about 25 K; the solver and the network agree to a few
        # 1e-3 K, the remainder being the diversion crossflow over the hydrodynamic development
        # length, which the network does not model. 0.004 K is less than a thousandth of the
        # spread, so a doubled or halved gap conductance (about 4 K) fails
        self.assertGreater(spread, 20.0)
        self.assertLess(
            error, 0.004, f"Largest difference from the network: {error:.3e} K"
        )

if __name__ == "__main__":
    unittest.main(verbosity=2)

(moose/modules/subchannel/verification/radial_conduction_verification/radial_conduction_7pin.i)

# Radial heat conduction verification in a 7-pin triangular assembly. The pins are unheated and
# bare, turbulent mixing is disabled (beta = 0, CT = 0), and a large lateral form loss
# coefficient Kij suppresses diversion crossflow. The duct heats the subchannels on one side of
# the assembly and cools those on the other side with the same total power, so heat crosses the
# assembly only by radial conduction through the gaps. The axial flow is low so that the
# temperature differences between subchannels are large. The outlet subchannel temperatures are
# compared with an independent solution of the same gap conductance network in
# test_radial_conduction.py.
T_in = 660 # K
mass_flux_in = 5.0 # kg/m2/s
P_out = 2.0e5 # Pa
duct_flux = 2.0e4 # W/m2

[GlobalParams]
  nrings = 2
  n_cells = 20
  flat_to_flat = 0.035
  heated_length = 0.5
  pitch = 0.012
[]

[TriSubChannelMesh]
  [subchannel]
    type = SCMTriAssemblyMeshGenerator
    pin_diameter = 0.01
    dwire = 0.0
    hwire = 0.0
    spacer_z = '0'
    spacer_k = '0'
    Kij = 1e8
  []
  [duct]
    type = SCMTriDuctMeshGenerator
    input = subchannel
  []
[]

# Constant properties close to sodium at 660 K, so that the network solution has constant
# conductances. With cp = cv, zero thermal expansion and porepressure_coefficient = 0,
# h = cp * T exactly
[FluidProperties]
  [sodium]
    type = SimpleFluidProperties
    density0 = 850
    thermal_expansion = 0
    cp = 1270
    cv = 1270
    thermal_conductivity = 70
    viscosity = 2.5e-4
    porepressure_coefficient = 0
  []
[]

[SubChannel]
  type = TriSubChannel1PhaseProblem
  fp = sodium
  P_out = ${P_out}
  compute_density = false
  compute_viscosity = false
  compute_power = true
  implicit = true
  segregated = true
  rtol = 1e-10
  atol = 1e-10
  friction_closure = 'Chen'
  mixing_closure = 'no_mixing'
  # The pin heat transfer closure only sets the pin surface temperature after the solve; with
  # q_prime = 0 the pins add no heat to the coolant. It is required whenever a pin mesh exists
  pin_HTC_closure = 'gnielinski'
  duct_HTC_closure = 'gnielinski'
[]

[SCMClosures]
  [Chen]
    type = SCMFrictionChenTodreas
  []
  [gnielinski]
    type = SCMHTCGnielinski
  []
  [no_mixing]
    type = SCMMixingConstantBeta
    beta = 0.0
    CT = 0.0
  []
[]

[Functions]
  # Heat the duct where x > 0 and cool it where x < 0. The duct nodes on x = 0 get no flux, so
  # the heated and cooled sides receive the same power
  [duct_flux_fn]
    type = ParsedFunction
    expression = 'if(x > 1e-6, ${duct_flux}, if(x < -1e-6, -${duct_flux}, 0))'
  []
[]

[ICs]
  [T_ic]
    type = ConstantIC
    variable = T
    value = ${T_in}
  []
  [P_ic]
    type = ConstantIC
    variable = P
    value = 0.0
  []
  [mu_ic]
    type = ViscosityIC
    variable = mu
    p = ${P_out}
    T = T
    fp = sodium
  []
  [rho_ic]
    type = RhoFromPressureTemperatureIC
    variable = rho
    p = ${P_out}
    T = T
    fp = sodium
  []
  [h_ic]
    type = SpecificEnthalpyFromPressureTemperatureIC
    variable = h
    p = ${P_out}
    T = T
    fp = sodium
  []
  [mdot_ic]
    type = ConstantIC
    variable = mdot
    value = 0.0
  []
  [q_prime_ic]
    type = ConstantIC
    variable = q_prime
    value = 0.0
  []
  [Dpin_ic]
    type = ConstantIC
    variable = Dpin
    value = 0.01
  []
  [duct_heat_flux_ic]
    type = FunctionIC
    variable = duct_heat_flux
    function = duct_flux_fn
  []
[]

[AuxKernels]
  [T_in_bc]
    type = ConstantAux
    variable = T
    boundary = inlet
    value = ${T_in}
    execute_on = 'timestep_begin'
  []
  [mdot_in_bc]
    type = SCMMassFlowRateAux
    variable = mdot
    boundary = inlet
    area = S
    mass_flux = ${mass_flux_in}
    execute_on = 'timestep_begin'
  []
[]

[VectorPostprocessors]
  # Sampled at every axial level for the independent network solution in
  # test_radial_conduction.py
  [channels]
    type = NodalValueSampler
    variable = 'T mdot S'
    block = subchannel
    sort_by = id
    execute_on = timestep_end
  []
  [duct]
    type = NodalValueSampler
    variable = 'duct_heat_flux'
    block = duct
    sort_by = id
    execute_on = timestep_end
  []
[]

[Postprocessors]
  [Total_power_SCMTHPowerPostprocessor]
    type = SCMTHPowerPostprocessor
  []
  [Total_power_SCMDuctPowerPostprocessor]
    type = SCMDuctHeatRatePostprocessor
  []
[]

[Executioner]
  type = Steady
[]

[Outputs]
  csv = true
[]

(moose/modules/subchannel/verification/radial_conduction_verification/test_radial_conduction.py)

"""
Radial heat conduction verification for the triangular-lattice subchannel solver.

radial_conduction_7pin.i heats the duct on one side of a 7-pin assembly and cools it on the other
with the same power. The pins are unheated, turbulent mixing is off, the fluid properties are
constant, and h = cp * T, so heat crosses the assembly only through the gap conductances.

This test solves the same conductance network independently of the solver and compares the
outlet subchannel temperatures. The network uses the documented radial conduction model of
TriSubChannel1PhaseProblem: each gap between channels i and j carries
k * w_ij * F_ij * dz / d_ij * (T_i - T_j), with F_ij = 0.66 (p / D) (w_ij / D)^-0.3. The gap width
w_ij is p - D for gaps between pins and the duct-to-pin gap for gaps facing the duct, and the
centroid distance d_ij is p / sqrt(3) for gaps involving an interior channel and p otherwise. The
duct adds q'' times the peripheral duct width of each edge or corner channel. The network is
marched in z with the solver's face mass flow rates, implicitly in the conduction term, which is
the discrete model the solver assembles. Only the channel positions, flow areas and mass flow
rates are taken from the solver output; the gap connectivity is rebuilt from the lattice.
"""

import csv
import math
import os
import subprocess
import unittest

import numpy as np


class TestRadialConductionVerification(unittest.TestCase):

    # Problem parameters (must match radial_conduction_7pin.i)
    T_in = 660.0
    cp = 1270.0
    k = 70.0
    pitch = 0.012
    pin_diameter = 0.01
    flat_to_flat = 0.035
    nrings = 2

    @classmethod
    def setUpClass(cls):
        cls.test_dir = os.path.dirname(os.path.abspath(__file__))
        cls.app = os.environ.get(
            "MOOSE_PYTHONUNITTEST_EXECUTABLE",
            os.path.join(cls.test_dir, "..", "..", "subchannel-opt"),
        )
        cls.input = os.path.join(cls.test_dir, "radial_conduction_7pin.i")
        cls.duct_to_pin_gap = 0.5 * (
            cls.flat_to_flat
            - (cls.nrings - 1) * cls.pitch * math.sqrt(3.0)
            - cls.pin_diameter
        )

    def _run(self, outbase, extra_args=None):
        """Run the solver and return the channel and duct samples at every axial level."""
        cmd = [
            self.app,
            "-i",
            self.input,
            f"Outputs/file_base={outbase}",
            "--no-color",
        ] + (extra_args or [])
        proc = subprocess.run(cmd, capture_output=True, text=True, cwd=self.test_dir)
        if proc.returncode != 0:
            self.fail(
                f"Solver failed ({outbase}):\n"
                f"--- stdout ---\n{proc.stdout[-3000:]}\n"
                f"--- stderr ---\n{proc.stderr[-1000:]}"
            )
        with open(os.path.join(self.test_dir, f"{outbase}_channels_0001.csv")) as fh:
            channels = list(csv.DictReader(fh))
        with open(os.path.join(self.test_dir, f"{outbase}_duct_0001.csv")) as fh:
            duct = list(csv.DictReader(fh))
        return channels, duct

    def _network_outlet_temperature(self, channels, duct):
        """Solve the gap conductance network; return (positions, T_solver, T_network)."""
        p, D, w_duct = self.pitch, self.pin_diameter, self.duct_to_pin_gap

        nodes = {}
        for r in channels:
            key = (round(float(r["x"]), 6), round(float(r["y"]), 6))
            nodes.setdefault(key, {})[round(float(r["z"]), 9)] = r
        positions = sorted(nodes)
        zs = sorted(nodes[positions[0]])
        n = len(positions)
        xy = np.array(positions)

        # Interior channels lie inside the outer pin ring; of the peripheral channels, the corner
        # channels have the smallest flow area
        area = np.array([float(nodes[c][zs[0]]["S"]) for c in positions])
        radius = np.hypot(xy[:, 0], xy[:, 1])
        r_outer_pins = (self.nrings - 1) * p
        kind = [
            (
                "interior"
                if r < r_outer_pins
                else ("corner" if math.isclose(s, area.min(), rel_tol=1e-9) else "edge")
            )
            for r, s in zip(radius, area)
        ]

        # Two channels share a gap when their centroids are the closest pair of their kinds.
        # Interior channels connect to interior and edge channels, edge channels to corner
        # channels; in a 7-pin assembly no two edge channels share a gap
        allowed = {("interior", "interior"), ("edge", "interior"), ("corner", "edge")}
        pairs = {}
        for i in range(n):
            for j in range(i + 1, n):
                kinds = tuple(sorted((kind[i], kind[j])))
                if kinds in allowed:
                    pairs.setdefault(kinds, []).append((math.dist(xy[i], xy[j]), i, j))
        gaps = []
        for kinds, candidates in pairs.items():
            d_min = min(c[0] for c in candidates)
            gaps += [(i, j, kinds) for d, i, j in candidates if d < 1.05 * d_min]

        # Conductance per unit length of each gap
        G = np.zeros((n, n))
        for i, j, kinds in gaps:
            if "interior" in kinds:
                width, dist = p - D, p / math.sqrt(3.0)
            else:
                width, dist = w_duct, p
            shape_factor = 0.66 * (p / D) * (width / D) ** -0.3
            g = self.k * width * shape_factor / dist
            G[i, i] += g
            G[j, j] += g
            G[i, j] -= g
            G[j, i] -= g

        # Duct heat per unit length: flux at the nearest duct node times the peripheral width
        duct_by_z = {}
        for r in duct:
            duct_by_z.setdefault(round(float(r["z"]), 9), []).append(r)
        corner_width = 2.0 / math.sqrt(3.0) * (D / 2.0 + w_duct)

        def duct_heat(i, z):
            if kind[i] == "interior":
                return 0.0
            x, y = positions[i]
            nearest = min(
                duct_by_z[z],
                key=lambda r: math.hypot(float(r["x"]) - x, float(r["y"]) - y),
            )
            width = p if kind[i] == "edge" else corner_width
            return float(nearest["duct_heat_flux"]) * width

        def mdot(z):
            return np.array([float(nodes[c][z]["mdot"]) for c in positions])

        h = np.full(n, self.cp * self.T_in)
        for iz in range(1, len(zs)):
            dz = zs[iz] - zs[iz - 1]
            q = np.array(
                [
                    0.5 * (duct_heat(i, zs[iz - 1]) + duct_heat(i, zs[iz])) * dz
                    for i in range(n)
                ]
            )
            h = np.linalg.solve(
                np.diag(mdot(zs[iz])) + G * dz / self.cp, mdot(zs[iz - 1]) * h + q
            )

        T_solver = np.array([float(nodes[c][zs[-1]]["T"]) for c in positions])
        return positions, T_solver, h / self.cp

    def test_matches_conductance_network(self):
        """Outlet temperatures agree with the independent network solution."""
        positions, T_solver, T_network = self._network_outlet_temperature(
            *self._run("radial_conduction_7pin_out")
        )
        spread = np.abs(T_solver - self.T_in).max()
        error = np.abs(T_solver - T_network).max()

        print(f"\nRadial conduction verification, 7-pin assembly")
        print(f'{"x [m]":>9} {"y [m]":>9} {"T solver [K]":>13} {"T network [K]":>14}')
        for (x, y), a, b in sorted(zip(positions, T_solver, T_network)):
            print(f"{x:+9.5f} {y:+9.5f} {a:13.5f} {b:14.5f}")
        print(f"Largest |T - T_in|: {spread:.3f} K, largest difference: {error:.2e} K")

        # Write the comparison for the verification page
        with open(
            os.path.join(self.test_dir, "radial_conduction_7pin_comparison.csv"), "w"
        ) as fh:
            fh.write("x,y,T_solver,T_network\n")
            for (x, y), a, b in sorted(zip(positions, T_solver, T_network)):
                fh.write(f"{x:.6f},{y:.6f},{a:.6f},{b:.6f}\n")

        # The heated and cooled sides span about 25 K; the solver and the network agree to a few
        # 1e-3 K, the remainder being the diversion crossflow over the hydrodynamic development
        # length, which the network does not model. 0.004 K is less than a thousandth of the
        # spread, so a doubled or halved gap conductance (about 4 K) fails
        self.assertGreater(spread, 20.0)
        self.assertLess(
            error, 0.004, f"Largest difference from the network: {error:.3e} K"
        )


if __name__ == "__main__":
    unittest.main(verbosity=2)

Test Description
Results
Input file