Conjugate Heat Transfer for Turbulent Channel Flow

In this tutorial, you will learn how to:

  • Couple NekRS to MOOSE for Conjugate Heat Transfer (CHT) for turbulent flow in a Tri-Structural Isotropic (TRISO) fuel unit cell

  • Couple NekRS to MOOSE for solving NekRS in non-dimensional form

  • Extract an initial condition from a NekRS restart file

  • Swap out NekRS for a different thermal-fluid MOOSE application, the Thermal Hydraulics Module (THM)

  • Compute heat transfer coefficients with NekRS

To access this tutorial,


cd cardinal/tutorials/gas_compact_cht

This tutorial also requires you to download some mesh files from Box. Please download the files from the gas_compact_cht folder here and place these files within the same directory structure in tutorials/gas_compact_cht.

This tutorial was developed with support from the NEAMS Thermal Fluids Center of Excellence and the NRIC Virtual Test Bed (VTB). You can find additional context on this model in our journal article Novak et al. (2022).

commentnote:Computing Needs

This tutorial requires HPC resources to run the NekRS cases. Please still read this tutorial if you do not have resources, because you will be able to run the THM cases and much of the Cardinal usage is agnostic of the particular thermal code being used.

Geometry and Computational Model

The geometry consists of a unit cell of a TRISO-fueled gas reactor compact INL (2016). A top-down view of the geometry is shown in Figure 1. The fuel is cooled by helium flowing in a cylindrical channel of diameter . Cylindrical fuel compacts containing randomly-dispersed TRISO particles at 15% packing fraction are arranged around the coolant channel in a triangular lattice. The TRISO particles use a conventional design that consists of a central fissile uranium oxycarbide kernel enclosed in a carbon buffer, an inner Pyrolitic Carbon (PyC) layer, a silicon carbide layer, and finally an outer PyC layer. The geometric specifications are summarized in Table 1. Heat is produced in the TRISO particles to yield a total power of 38 kW.

Figure 1: TRISO-fueled gas reactor compact unit cell

Table 1: Geometric specifications for a TRISO-fueled gas reactor compact

ParameterValue (cm)
Coolant channel diameter, 1.6
Fuel compact diameter, 1.27
Fuel-to-coolant center distance, 1.628
Height160
TRISO kernel radius214.85e-4
Buffer layer radius314.85e-4
Inner PyC layer radius354.85e-4
Silicon carbide layer radius389.85e-4
Outer PyC layer radius429.85e-4

Heat Conduction Model

The MOOSE heat transfer module is used to solve for energy conservation in the solid, with the time derivative neglected in order to more quickly approach steady state. The solid mesh is shown in Figure 2. The TRISO particles are homogenized into the compact regions - all material properties in the heterogeneous regions are taken as volume averages of the various constituent materials.

Figure 2: Mesh for the solid heat conduction model

The volumetric power density is set to a sinusoidal function of ,

where is a coefficient used to ensure that the total power produced is 38 kW and is the domain height. The power is uniform in the radial direction in each compact. On the coolant channel surface, a Dirichlet temperature is provided by NekRS/THM. All other boundaries are insulated. We will run the solid model first, so we must specify an initial condition for the wall temperature, which we simply set to a linear variation between the inlet and the outlet based on the nominal temperature rise.

NekRS Model

NekRS solves the incompressible k-tau RANS equations. To relate to , is selected. The resolution of the mesh near all no-slip boundaries ensures that such that the NekRS model is a wall-resolved model.

Here, the NekRS case will be set up in non-dimensional form. This is a convenient technique for fluid solvers, because the solution will be correct for any flow at a given Reynolds/Peclet/etc. number in this geometry. A common user operation is also to "ramp" up a CFD solve for stability purposes. Solving in non-dimensional form means that you do not need to also "ramp" absolute values for velocity/pressure/temperature/etc. - switching to a different set of non-dimensional numbers is as simple as just modifying the thermophysical "properties" of the fluid (e.g. viscosity, conductivity, etc.).

When the Navier-Stokes equations are written in non-dimensional form, the "density" in the [VELOCITY] block becomes unity because

(1)

The "viscosity" becomes the coefficient on the viscous stress term (, where is the Reynolds number). In NekRS, specifying diffusivity = -50.0 is equivalent to specifying diffusivity = 0.02 (i.e. ), or a Reynolds number of 50.0.

In non-dimensional form, the rhoCp term in the [TEMPERATURE] block becomes unity because

(2)

The "conductivity" indicates the coefficient on the diffusion kernel, which in non-dimensional form is equal to , where is the Peclet number. In NekRS, specifying conductivity = -35 is equivalent to specifying conductivity = 0.02857 (i.e. ), or a Peclet number of 35.

The inlet mass flowrate is 0.0905 kg/s; with the channel diameter of 1.6 cm and material properties of helium, this results in a Reynolds number of 223214 and a Prandtl number of 0.655. This highly-turbulent flow results in extremely thin momentum and thermal boundary layers on the no-slip surfaces forming the periphery of the coolant channel. In order to resolve the near-wall behavior, an extremely fine mesh is required in the NekRS simulation. To accelerate the overall coupled CHT case that is of interest in this tutorial, the NekRS model is split into a series of calculations:

  1. We first run a partial-height, periodic flow-only case to obtain converged , , and distributions.

  2. Then, we extrapolate the and to the full-height case.

  3. Finally, we use the converged, full-height and distributions to transport a temperature passive scalar in a CHT calculation with MOOSE.

As a rough estimate, solving the coupled mass-momentum equations requires about an order of magnitude more compute time than a passive scalar equation in NekRS. Therefore, by "freezing" the and solutions, we can dramatically reduce the total cost of the coupled CHT calculation.

The periodic flow model has a height of and takes about 10,000 core hours to reach steady state. The extrapolation from to then is a simple postprocessing operation. Because steps 1 and 2 are done exclusively using NekRS, and the periodic flow solve takes a very long time to run (from the perspective of a tutorial, at least), we omit steps 1 and 2 from this tutorial. Instead, we begin straight away from a full-height periodic restart file with the - RANS model. This restart file is the periodic_flow_solve.f00001 file. If you open the restart file, you will see the following converged NekRS predictions of velocity, , and . Figure 4 shows the axial velocity, , and on a line plot through the center of the channel.

Figure 3: Converged NekRS predictions of velocity magnitude (left), (middle), and (right) for the periodic case (all in non-dimensional form)

Figure 4: Converged NekRS predictions of velocity magnitude (left), (middle), and (right) for the periodic case along the channel mide-plane (all in non-dimensional form)

For the conjugate heat transfer case, we will load this restart file, compute from the loaded solutions for and , and then transport temperature with coupling to MOOSE heat conduction. The NekRS files are:

  • ranstube.re2: NekRS mesh

  • ranstube.par: High-level settings for the solver, boundary condition mappings to sidesets, and the equations to solve

  • ranstube.udf: User-defined C++ functions for on-line postprocessing and model setup

  • ranstube.oudf: User-defined OCCA kernels for boundary conditions and source terms

The NekRS mesh is shown in Figure 5. Boundary 1 is the inlet, boundary 2 is the outlet, and boundary 3 is the wall. The same mesh was used for the periodic flow solve, except with a shorter height.

Figure 5: Mesh for the NekRS RANS model

Next, the .par file contains problem setup information. This input sets up a nondimensional passive scalar solution, loading , , , and from a restart file. In order to "freeze," or turn off the , , , and solves, we set solver = none in the [VELOCITY], [SCALAR01] ( passive scalar), and [SCALAR02] ( passive scalar) blocks. The only equation that NekRS will solve is for temperature.

# Model of a turbulent channel flow.
#
# L_ref = 0.016 m          coolant channel diameter
# T_ref = 598 K            coolant inlet temperature
# dT_ref = 82.89 K         coolant nominal temperature rise in the unit cell
# U_ref = 81.089           coolant inlet velocity
# length = 1.6 m           coolant channel flow length
#
# rho_0 = 5.5508 kg/m3     coolant density
# mu_0 = 3.22639e-5 Pa-s   coolant dynamic viscosity
# Cp_0 = 5189 J/kg/K       coolant isobaric heat capacity
# k_0 = 0.2556 W/m/K       coolant thermal conductivity
#
# Re = 223214
# Pr = 0.655
# Pe = 146205

[GENERAL]
  startFrom = periodic_flow_solve.f00001

  polynomialOrder = 7
  dt = 6e-3
  timeStepper = BDF2

  stopAt = numSteps
  numSteps = 1000

  # write an output file every specified time steps
  writeControl = steps
  writeInterval = 1000

[PROBLEMTYPE]
  stressFormulation = true

[PRESSURE]
  residualTol = 1e-04

[VELOCITY]
  solver = none
  viscosity = -223214
  density = 1.0
  boundaryTypeMap = v, o, W
  residualTol = 1e-06

[TEMPERATURE]
  rhoCp = 1.0
  conductivity = -146205
  boundaryTypeMap = t, o, f
  residualTol = 1e-06

[SCALAR01] # k
  solver = none
  boundaryTypeMap = t, I, o
  residualTol = 1e-06
  rho = 1.0
  diffusivity = -223214

[SCALAR02] # tau
  solver = none
  boundaryTypeMap = t, I, o
  residualTol = 1e-06
  rho = 1.0
  diffusivity = -223214
(tutorials/gas_compact_cht/ranstube.par)

Next, the .udf file is used to setup initial conditions and define how should be computed based on and the restart values of and . In turbulent_props, a user-defined function, we use from the input file in combination with the and (read from the restart file later in the .udf file) to adjust the total diffusion coefficient on temperature to according to the definition of turbulent Prandtl number. This adjustment must happen on device, in a new GPU kernel we name scalarScaledAddKernel. This kernel will be defined in the .oudf file; we instruct the JIT compilation to compile this new kernel by calling udfBuildKernel.

Then, in UDF_Setup we set an initial condition for fluid temperature (the first scalar in the nrs->cds->S array that holds all the scalars). In this function, we also store the value of computed in the restart file.

#include <math.h>
#include "udf.hpp"

occa::memory o_mu_t_from_restart;

#define length 100.0    // non-dimensional channel length
#define Pr_T 0.91       // turbulent Prandtl number

static occa::kernel scalarScaledAddKernel;

void turbulent_props(nrs_t *nrs, double time, occa::memory o_U, occa::memory o_S,
            occa::memory o_UProp, occa::memory o_SProp)
{
  // fetch the laminar thermal conductivity and specify the desired turbulent Pr
  dfloat k_laminar;
  platform->options.getArgs("SCALAR00 DIFFUSIVITY", k_laminar);

  // o_diff is the turbulent conductivity for all passive scalars, so we grab the first
  // slice, since temperature is the first passive scalar
  occa::memory o_mu = nrs->cds->o_diff + 0 * nrs->cds->fieldOffset[0] * sizeof(dfloat);

  scalarScaledAddKernel(0, nrs->cds->mesh[0]->Nlocal, k_laminar, 1.0 / Pr_T,
    o_mu_t_from_restart /* turbulent viscosity */, o_mu /* laminar viscosity */);
}

void UDF_LoadKernels(occa::properties & kernelInfo)
{
  scalarScaledAddKernel = oudfBuildKernel(kernelInfo, "scalarScaledAdd");
}

void UDF_Setup(nrs_t *nrs)
{
  mesh_t *mesh = nrs->cds->mesh[0];
  int n_gll_points = mesh->Np * mesh->Nelements;

  // allocate space for the k and tau that we will read from the file
  dfloat * mu_t = (dfloat *) calloc(n_gll_points, sizeof(dfloat));

  for (int n = 0; n < n_gll_points; ++n)
  {
    // initialize the temperature to a linear variation from 0 to 1
    dfloat z = mesh->z[n];
    nrs->cds->S[n + 0 * nrs->cds->fieldOffset[0]] = z / length;

    // fetch the restart value of mu_T, which is equal to k * tau, the second and third
    // passive scalars
    mu_t[n] = nrs->cds->S[n + 1 * nrs->cds->fieldOffset[0]] *
              nrs->cds->S[n + 2 * nrs->cds->fieldOffset[0]];
  }

  o_mu_t_from_restart = platform->device.malloc(n_gll_points * sizeof(dfloat), mu_t);

  udf.properties = &turbulent_props;
}

void UDF_ExecuteStep(nrs_t * nrs, double time, int tstep)
{
}
(tutorials/gas_compact_cht/ranstube.udf)

In the .oudf file, we define boundary conditions for temperature and also the form of the scalarScaledAdd kernel that we use to compute . The inlet boundary is set to a temperature of 0 (a dimensional temperature of ), while the fluid-solid interface will receive a heat flux from MOOSE.

@kernel void scalarScaledAdd(const dlong offset, const dlong N,
                             const dfloat a,
                             const dfloat b,
                             @restrict const dfloat* X,
                             @restrict dfloat* Y)
{
  for (dlong n = 0; n < N; ++n; @tile(256,@outer,@inner))
    if (n < N)
      Y[n] = a + b * X[n + offset];
}

// Boundary conditions are only needed for temperature, since the solves of pressure,
// velocity, k, and tau are all frozen

void scalarDirichletConditions(bcData *bc)
{
  // inlet temperature
  bc->s = 0.0;
}

void scalarNeumannConditions(bcData *bc)
{
  // wall (temperature)
  // note: when running with Cardinal, Cardinal will allocate the usrwrk
  // array. If running with NekRS standalone (e.g. nrsmpi), you need to
  // replace the usrwrk with some other value or allocate it youself from
  // the udf and populate it with values.
  bc->flux = bc->usrwrk[bc->idM];
}
(tutorials/gas_compact_cht/ranstube.oudf)

THM Model

THM solves the 1-D area averages of the Navier-Stokes equations. The Churchill correlation is used for and the Dittus-Boelter correlation is used for Berry et al. (2014).

The THM mesh contains 150 elements; the mesh is constucted automatically within THM. The heat flux imposed in the THM elements is obtained by area averaging the heat flux from the heat conduction model in 150 layers along the fluid-solid interface. For the reverse transfer, the wall temperature sent to MOOSE heat conduction is set to a uniform value along the fluid-solid interface according to a nearest-node mapping to the THM elements.

Because the thermal-fluids tools run after the MOOSE heat conduction application, initial conditions are only required for pressure, fluid temperature, and velocity, which are set to uniform distributions.

CHT Coupling

In this section, MOOSE, NekRS, and THM are coupled for CHT. Two separate simulations are performed here:

  • Coupling of NekRS with MOOSE

  • Coupling of THM with MOOSE

By individually describing the two setups, you will understand that NekRS is really as interchangeable as any other MOOSE application for coupling. In addition, we will describe how to generate a heat transfer coefficient from the NekRS-MOOSE coupling to input into the THM-MOOSE coupling to demonstrate an important multiscale use case of NekRS.

NekRS-MOOSE

In this section, we describe the coupling of MOOSE and NekRS.

Solid Input Files

The solid phase is solved with the MOOSE heat transfer module, and is described in the solid_nek.i input. We define a number of constants at the beginning of the file and set up the mesh.

# This input models conjugate heat transfer between MOOSE and NekRS for
# a single coolant channel. This input file should be run with:
#
# cardinal-opt -i common_input.i solid_nek.i

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}
unit_cell_mdot = ${fparse mdot / (n_bundles * n_coolant_channels_per_block)}

# compute the volume fraction of each TRISO layer in a TRISO particle
# for use in computing average thermophysical properties
kernel_fraction = ${fparse kernel_radius^3 / oPyC_radius^3}
buffer_fraction = ${fparse (buffer_radius^3 - kernel_radius^3) / oPyC_radius^3}
ipyc_fraction = ${fparse (iPyC_radius^3 - buffer_radius^3) / oPyC_radius^3}
sic_fraction = ${fparse (SiC_radius^3 - iPyC_radius^3) / oPyC_radius^3}
opyc_fraction = ${fparse (oPyC_radius^3 - SiC_radius^3) / oPyC_radius^3}

# multiplicative factor on assumed heat source distribution to get correct magnitude
q0 = ${fparse unit_cell_power / (4.0 * unit_cell_height * compact_diameter * compact_diameter / 4.0)}

# Time step interval on which to exchange data between NekRS and MOOSE
N = 50
U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3

[Mesh]
  type = FileMesh
  file = solid_rotated.e
[]
(tutorials/gas_compact_cht/solid_nek.i)

Next, we define a temperature variable T, and specify the governing equations and boundary conditions we will apply.

[Variables]
  [T]
    initial_condition = ${inlet_T}
  []
[]

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = 'compacts'
  []
[]

[BCs]
  [pin_outer]
    type = MatchedValueBC
    variable = T
    v = fluid_temp
    boundary = 'fluid_solid_interface'
  []
[]

[ICs]
  [power]
    type = FunctionIC
    variable = power
    function = axial_power
    block = 'compacts'
  []
  [fluid_temp]
    type = FunctionIC
    variable = fluid_temp
    function = axial_fluid_temp
  []
[]
(tutorials/gas_compact_cht/solid_nek.i)

The MOOSE heat transfer module will receive a wall temperature from NekRS in the form of an AuxVariable, so we define a receiver variable for the temperature, as fluid_temp. The MOOSE heat transfer module will also send heat flux to NekRS, which we compute as another AuxVariable named flux, which we compute with a DiffusionFluxAux auxiliary kernel. Finally, we define another auxiliary variable for the imposed power, which we will not receive from a coupled application, but instead set within the solid input file.

[AuxVariables]
  [fluid_temp]
  []
  [flux]
    family = MONOMIAL
    order = CONSTANT
  []
  [power]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [flux]
    type = DiffusionFluxAux
    diffusion_variable = T
    component = normal
    diffusivity = thermal_conductivity
    variable = flux
    boundary = 'fluid_solid_interface'
  []
[]
(tutorials/gas_compact_cht/solid_nek.i)

Next, we use functions to define the thermal conductivities. The material properties for the TRISO compacts are taken as volume averages of the various constituent materials. We also use functions to set initial conditions for power (a sinusoidal function) and the initial fluid temperature (a linear variation from inlet to outlet).

[Functions]
  [k_graphite]
    type = ParsedFunction
    expression = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction
    expression = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction
    expression = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names = 'k_TRISO k_graphite'
    symbol_values = 'k_TRISO k_graphite'
  []
  [axial_power] # volumetric power density
    type = ParsedFunction
    expression = 'sin(pi * z / ${unit_cell_height}) * ${q0}'
  []
  [axial_fluid_temp]
    type = ParsedFunction
    expression = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / ${unit_cell_mdot} / ${fluid_Cp}'
  []
[]

[Materials]
  [graphite]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_graphite
    temp = T
    block = 'graphite'
  []
  [compacts]
    type = HeatConductionMaterial
    thermal_conductivity_temperature_function = k_compacts
    temp = T
    block = 'compacts'
  []
[]
(tutorials/gas_compact_cht/solid_nek.i)

We define a number of postprocessors for querying the solution as well as for normalizing the heat flux.

[Postprocessors]
  [flux_integral]
    # evaluate the total heat flux for normalization
    type = SideDiffusiveFluxIntegral
    diffusivity = thermal_conductivity
    variable = T
    boundary = 'fluid_solid_interface'
  []
  [max_fuel_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'compacts'
  []
  [max_block_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = 'graphite'
  []
  [max_flux]
    type = ElementExtremeValue
    variable = flux
  []
  [synchronization_to_nek]
    type = Receiver
    default = 1.0
  []

  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = 'compacts'
  []
[]
(tutorials/gas_compact_cht/solid_nek.i)

For visualization purposes only, we add LayeredAverages to average the fuel and block temperatures in layers in the direction. We also add a LayeredSideAverage to average the heat flux along a boundary, again in layers in the direction. We then output the results of these userobjects to CSV using SpatialUserObjectVectorPostprocessors and by setting csv = true in the output.

[UserObjects]
  [average_fuel_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_block_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'graphite'
  []
  [average_flux_axial]
    type = LayeredSideAverage
    variable = flux
    direction = z
    num_layers = ${num_layers_for_plots}
    boundary = 'fluid_solid_interface'
  []
[]

[VectorPostprocessors]
  [fuel_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fuel_axial
  []
  [block_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_block_axial
  []
  [flux_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_flux_axial
  []
[]

[Outputs]
  exodus = true
  print_linear_residuals = false
  hide = 'synchronization_to_nek'

  [csv]
    file_base = 'csv/solid_nek'
    type = CSV
  []
[]
(tutorials/gas_compact_cht/solid_nek.i)

Finally, we add a TransientMultiApp that will run a MOOSE-wrapped NekRS simulation. Then, we add four different transfers to/from NekRS:

[MultiApps]
  [nek]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'nek.i'
    sub_cycling = true
    execute_on = timestep_end
  []
[]

[Transfers]
  [heat_flux_to_nek]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = flux
    variable = avg_flux
    from_boundaries = 'fluid_solid_interface'
    to_boundaries = '3'
    to_multi_app = nek
  []
  [flux_integral_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = flux_integral
    from_postprocessor = flux_integral
    to_multi_app = nek
  []
  [temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    variable = fluid_temp
    from_boundaries = '3'
    to_boundaries = 'fluid_solid_interface'
    from_multi_app = nek
  []
  [synchronization_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = transfer_in
    from_postprocessor = synchronization_to_nek
    to_multi_app = nek
  []
[]
(tutorials/gas_compact_cht/solid_nek.i)

Finally, although our heat conduction model does not have any time derivatives (in order to more quickly reach steady state), NekRS must be run as a transient. In order to control the interval on which NekRS couples to MOOSE, we use a MOOSE time step that is N times bigger than the NekRS time step, taking care to account for the fact that the time step in the ranstube.par file is actually a nondimensional time step. By setting subcycling = true for the TransientMultiApp, MOOSE is only run every N time steps (the subcycling setting), and that the NekRS solution is only copied to/from GPU every N time steps (the minimize transfers setting).

Finally, we will run the MOOSE heat conduction model until the relative change in the solution (across the entire domain) differs by less than from the previous time step. While this may seem like a large tolerance, because the difference is taken over the entire domain, you will see changes only on the order of the fourth or fifth decimal place in temperatures.

[Executioner]
  type = Transient
  dt = '${fparse N * nek_dt * t0}'

  nl_abs_tol = 1e-5
  nl_rel_tol = 1e-16
  petsc_options_value = 'hypre boomeramg'
  petsc_options_iname = '-pc_type -pc_hypre_type'

  steady_state_detection = true
  steady_state_tolerance = 1e-1
[]
(tutorials/gas_compact_cht/solid_nek.i)

Fluid Input Files

The Nek wrapping is described in the nek.i input file. We first define a few file-local variables, and then build a mesh mirror with a NekRSMesh. By setting boundary = '3', we indicate that boundary 3 will be coupled via CHT. In this example, we don't have any volume-based coupling, but we set volume = true to visualize the NekRS solution in Exodus format (as a first-order interpolation of the 7-th order polynomial solution).

When solving in MOOSE, all the coupled applications must either be in non-dimensional form, or all in dimensional form. In this tutorial, we have set up the MOOSE input file in dimensional form, but the NekRS case files were set up in non-dimensional form. In order for MOOSE's transfers to correctly find the closest nodes in the solid mesh to corresponding locations in NekRS, the entire mesh "mirror" through which data transfers occur must be scaled by a factor of to return to dimensional units (because the coupled MOOSE application is in dimensional units). This scaling is specified by the scaling parameter. In other words, the scaling indicates by what factor we should internally scale spatial lengths by to correctly send data between NekRS and MOOSE.

# copy-pasta from common_input.i

channel_diameter = 0.016                 # diameter of the coolant channels (m)
height = 6.343                           # height of the full core (m)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected
num_layers_for_plots = 50                # number of layers to average fields over for plotting

[Mesh]
  type = NekRSMesh
  boundary = '3'
  volume = true
  scaling = ${channel_diameter}
[]
(tutorials/gas_compact_cht/nek.i)

Next, we define additional parameters to describe how NekRS interacts with MOOSE with the NekRSProblem.

To allow conversion between a non-dimensional NekRS solve and a dimensional MOOSE coupled heat conduction application, the characteristic scales used to establish the non-dimensional problem are provided. Definitions for these non-dimensional scales are available here.

warningwarning

These characteristic scales are used by Cardinal to dimensionalize the NekRS solution into the units that the coupled MOOSE application expects. You still need to properly non-dimensionalize the NekRS input files (to be discussed later).

We also indicate that we are going to restrict the data copies to/from GPU to only occur on the time steps that NekRS is coupled to MOOSE.

[Problem]
  type = NekRSProblem
  casename = 'ranstube'
  has_heat_source = false

  nondimensional = true
  U_ref = '${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}'
  T_ref = ${inlet_T}
  dT_ref = '${fparse power / mdot / fluid_Cp * unit_cell_height / height}'
  L_ref = ${channel_diameter}
  rho_0 = ${fluid_density}
  Cp_0 = ${fluid_Cp}

  synchronization_interval = parent_app
[]
(tutorials/gas_compact_cht/nek.i)

For time stepping, we will allow NekRS to control its own time stepping by using the NekTimeStepper. We will output both Exodus and CSV format data as well.

[Executioner]
  type = Transient

  [TimeStepper]
    type = NekTimeStepper
    min_dt = 1e-10
  []
[]

[Outputs]
  exodus = true

  [screen]
    type = Console
    hide = 'boundary_flux inlet_T outlet_T max_T flux_integral transfer_in'
  []

  [csv]
    file_base = 'csv/nek'
    type = CSV
  []
[]
(tutorials/gas_compact_cht/nek.i)

Then, we define a few postprocessors to use for querying the solution. While we hide these from the console output above, they will still be output to CSV.

[Postprocessors]
  [boundary_flux]
    type = NekHeatFluxIntegral
    boundary = '3'
  []
  [inlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '1'
  []
  [outlet_T]
    type = NekSideAverage
    field = temperature
    boundary = '2'
  []
  [max_T]
    type = NekVolumeExtremeValue
    field = temperature
    value_type = max
  []
[]
(tutorials/gas_compact_cht/nek.i)

Finally, we add a few userobjects to compute the necessary terms in a heat transfer coefficient in a series of equally-spaced axial layers along the flow direction as

(3)

where is the layer index, is the average wall heat flux in layer , is the heat transfer coefficient in layer , is the average wall temperature in layer , and is the volume-averaged temperature in layer . All terms in Eq. (3) can be computed using objects in the MOOSE and NekRS-wrapped input files. In the NekRS input file, we add userobjects to compute average wall temperatures and bulk fluid temperatures in axial layers with NekBinnedSideAverage and NekBinnedVolumeAverage objects. We then output the results of these userobjects to CSV using SpatialUserObjectVectorPostprocessors. For the time being, we turn off the calculation of these objects until we are ready to describe the generation of heat transfer coefficients in Heat Transfer Coefficients .

[UserObjects]
  [average_fuel_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_block_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'graphite'
  []
  [average_flux_axial]
    type = LayeredSideAverage
    variable = flux
    direction = z
    num_layers = ${num_layers_for_plots}
    boundary = 'fluid_solid_interface'
  []
[]

[VectorPostprocessors]
  [fuel_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_fuel_axial
  []
  [block_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_block_axial
  []
  [flux_axial_avg]
    type = SpatialUserObjectVectorPostprocessor
    userobject = average_flux_axial
  []
[]
(tutorials/gas_compact_cht/solid_nek.i)

THM-MOOSE

In the previous section, we described a coupling of NekRS and MOOSE. In this section, we describe the coupling of MOOSE and THM, a lower-order fluid solver.

Solid Input Files

The solid phase is again solved with the MOOSE heat transfer module, and is described in the solid_thm.i input. The solid input file is almost exactly the same as the heat conduction model used for coupling to NekRS that has already been described, so we only focus on the differences. Instead of creating a NekRS sub-application, we create a THM sub-application.

[MultiApps]
  [thm]
    type = FullSolveMultiApp
    input_files = 'thm.i'
    execute_on = timestep_end
    max_procs_per_app = 1
    bounding_box_padding = '0.1 0.1 0'
  []
[]
(tutorials/gas_compact_cht/solid_thm.i)

Because THM is a 1-D code, we must perform an averaging operation on the surface heat flux computed by MOOSE before it is sent to THM. This requires us to use slightly different MOOSE transfers to couple THM and MOOSE. Instead of sending a 2-D surface flux distribution, we now compute the wall average heat flux in layers and send to THM with a MultiAppGeneralFieldUserObjectTransfer. To receive temperature from THM, we use the same MultiAppGeneralFieldNearestLocationTransfer that we used when coupling to NekRS - but now, we are coupling a 1-D model to a 3-D model instead of a 3-D model to a 3-D model. MOOSE's transfers are generally dimension agnostic.

[Transfers]
  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer
    variable = q_wall
    to_multi_app = thm
    source_user_object = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = T_wall
    from_multi_app = thm
    variable = fluid_temp
  []
[]
(tutorials/gas_compact_cht/solid_thm.i)

Other postprocessing features in this MOOSE model are largely the same as in the NekRS-MOOSE case; except to compute the wall heat flux to send to THM, we use a LayeredSideDiffusiveFluxAverage, which computes heat flux on a boundary given the temperature computed by MOOSE.

[UserObjects]
  [average_fuel_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'compacts'
  []
  [average_block_axial]
    type = LayeredAverage
    variable = T
    direction = z
    num_layers = ${num_layers_for_plots}
    block = 'graphite'
  []
  [q_wall_avg]
    type = LayeredSideDiffusiveFluxAverage
    boundary = fluid_solid_interface
    direction = z
    variable = T
    diffusivity = thermal_conductivity
    # Note: make this to match the num_elems in the channel
    num_layers = 50
  []
[]
(tutorials/gas_compact_cht/solid_thm.i)

Fluid Input Files

The fluid phase will be solved with THM, and is described in the thm.i file. The THM input file is built using syntax specific to THM - we will only briefly cover the syntax, and instead refer users to the THM manuals for more information. First, we define a number of constants at the beginning of the file and apply some global settings. We set the initial conditions for pressure, velocity, and temperature and indicate the fluid Equation of State (EOS) object using IdealGasFluidProperties.

# copy-pasta from common.i
inlet_T = 598.0                          # inlet fluid temperature (K)
mdot = ${fparse 117.3 / (12 * 108)}      # fluid mass flowrate (kg/s)
outlet_P = 7.1e6                         # fluid outlet pressure (Pa)
channel_diameter = 0.016                 # diameter of the coolant channels (m)
height = 1.6                             # height of the unit cell (m)
num_layers_for_plots = 50

[GlobalParams]
  initial_p = ${outlet_P}
  initial_T = ${inlet_T}
  initial_vel = 1.0

  rdg_slope_reconstruction = full
  closures = none
  fp = helium
[]

[Closures]
  [none]
    type = Closures1PhaseNone
  []
[]

[FluidProperties]
  [helium]
    type = IdealGasFluidProperties
    molar_mass = 4e-3
    gamma = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k = 0.2556
    mu = 3.22639e-5
  []
[]
(tutorials/gas_compact_cht/thm.i)

Next, we define the "components" in the domain. These components essentially consist of the physics equations and boundary conditions solved by THM, but expressed in THM-specific syntax. These components define single-phase flow in a pipe, an inlet mass flowrate boundary condition, an outlet pressure condition, and heat transfer to the pipe wall.

[Components]
  [channel]
    type = FlowChannel1Phase
    position = '0 0 0'
    orientation = '0 0 1'

    A = '${fparse pi * channel_diameter * channel_diameter / 4}'
    D_h = ${channel_diameter}
    length = ${height}
    n_elems = 50
  []

  [inlet]
    type = InletMassFlowRateTemperature1Phase
    input = 'channel:in'
    m_dot = ${mdot}
    T = ${inlet_T}
  []

  [outlet]
    type = Outlet1Phase
    input = 'channel:out'
    p = ${outlet_P}
  []

  [ht_ext]
    type = HeatTransferFromExternalAppHeatFlux1Phase
    flow_channel = channel
    P_hf = '${fparse channel_diameter * pi}'
  []
[]
(tutorials/gas_compact_cht/thm.i)

Associated with these components are a number of closures, defined as materials. We set up the Churchill correlation for the friction factor and the Dittus-Boelter correlation for the convective heat transfer coefficient. Additional materials are created to represent dimensionless numbers and other auxiliary terms, such as the wall temperature. As can be seen here, the Material system is not always used to represent quantities traditionally thought of as "material properties."

[Materials]
  # wall friction closure
  [f_mat]
    type = ADWallFrictionChurchillMaterial
    block = channel
    D_h = D_h
    f_D = f_D
    vel = vel
    rho = rho
    mu = mu
  []

  # Wall heat transfer closure (all important is in Nu_mat)
  [Re_mat]
    type = ADReynoldsNumberMaterial
    block = channel
    Re = Re
    D_h = D_h
    mu = mu
    vel = vel
    rho = rho
  []
  [Pr_mat]
    type = ADPrandtlNumberMaterial
    block = channel
    cp = cp
    mu = mu
    k = k
  []
  [Nu_mat]
    type = ADParsedMaterial
    block = channel
    # Dittus-Boelter
    expression = '0.023 * pow(Re, 0.8) * pow(Pr, 0.4)'
    property_name = 'Nu'
    material_property_names = 'Re Pr'
  []
  [Hw_mat]
    type = ADConvectiveHeatTransferCoefficientMaterial
    block = channel
    D_h = D_h
    Nu = Nu
    k = k
  []

  [T_wall]
    type = ADTemperatureWall3EqnMaterial
    Hw = Hw
    T = T
    q_wall = q_wall
  []
[]
(tutorials/gas_compact_cht/thm.i)

THM computes the wall temperature to apply a boundary condition in the MOOSE heat transfer module. To convert the T_wall material into an auxiliary variable, we use the ADMaterialRealAux.

[AuxVariables]
  [T_wall]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [Tw_aux]
    type = ADMaterialRealAux
    block = channel
    variable = T_wall
    property = T_wall
  []
[]
(tutorials/gas_compact_cht/thm.i)

Finally, we set the preconditioner, a Transient executioner, and set an Exodus output. We will run THM to convergence based on a tight steady state relative tolerance of .

[Preconditioning]
  [pc]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  dt = 0.1
  start_time = 0

  # end_time = 1000
  steady_state_detection = true
  steady_state_tolerance = 1e-08

  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-6

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  solve_type = NEWTON
  line_search = basic
[]

[Outputs]
  exodus = true
  print_linear_residuals = false

  [screen]
    type = Console
    outlier_variable_norms = false
    hide = 'max_p min_p'
  []

  [csv]
    file_base = 'csv/thm'
    type = CSV
  []
[]
(tutorials/gas_compact_cht/thm.i)

Execution and Postprocessing

To run the coupled NekRS-MOOSE calculation,


mpiexec -np 500 cardinal-opt -i common_input.i solid_nek.i

This will run with 500 MPI processes. To run the coupled THM-MOOSE calculation,


mpiexec -np 2 cardinal-opt -i common_input.i solid_thm.i --n-threads=2

which will run with 2 MPI ranks and 2 OpenMP threads per rank.

When the two simulations have completed, you will have created a number of different output files:

  • solid_nek_out.e, an Exodus file with the MOOSE solution for the NekRS-MOOSE calculation

  • solid_nek_out_nek0.e, an Exodus file with the NekRS solution that was ultimately transferred in/out of MOOSE

  • ranstube0.f<n>, a (custom format) NekRS output file with the NekRS solution; to convert to a format viewable in Paraview, consult the NekRS documentation

  • solid_thm_out.e, an Exodus file with the MOOSE solution for the THM-MOOSE calculation

  • solid_thm_out_thm0.e, an Exodus file with the THM solution

All CSV output files are placed in the csv directory. In discussing the results, we will refer to the nominal Dittus-Boelter THM simulations as the "baseline" THM-MOOSE simulations. In Heat Transfer Coefficients when we correct the Dittus-Boelter correlation based on the NekRS-MOOSE simulations, we will refer to those THM-MOOSE simulations are "corrected" THM-MOOSE simulations.

Figure 6 shows the solid temperature predicted for the NekRS-MOOSE simulation, the baseline THM-MOOSE simulation, and the difference between the two (NekRS minus THM). Temperatures are highest in the six compact regions, slightly downstream from the core midplane due to the combined effects of convective heat transfer and the imposed sinusoidal power distribution. The solid temperature distributions predicted with either NekRS fluid coupling or THM fluid coupling agree very well with each other - the maximum solid temperature difference is 16.3 K.

Figure 6: Solid temperature predicted for NekRS-MOOSE, baseline THM-MOOSE, and the difference between the two (NekRS minus baseline THM).

This code difference can be attributed to a number of different sources:

  • Because THM averages the solid heat flux along the coolant channel perimeter, THM is incapable of representing the azimuthal variation in heat flux that occurs due to the power production in the six fuel compacts.

  • NekRS uses a wall-resolved methodology, whereas THM uses the Dittus-Boelter correlation. Assuming the maximum temperature difference could be solely attributed to a difference in Nusselt number, a 16.3 K difference is equivalent to about a 10% difference in convective heat transfer coefficient. Errors as high as 25% may be observed when using a simple model such as the Dittus-Boelter correlation Incropera and Dewitt (2011), and a 10% difference is well within the realm of the data scatter used to originally fit the Dittus-Boelter model Winterton (1998)

Figure 7 shows the solid temperature predictions on the axial midplane for the NekRS-MOOSE, the baseline THM-MOOSE simulation, and the difference between the two (NekRS minus baseline THM). Figure 7 shows with greater clarity that the highest temperatures are observed in the fuel compacts, and that the inability of THM to resolve the angular variation in wall heat flux contributes to some spatial variation of the solid temperature difference along planes.

Figure 7: Solid temperature predicted along the axial mid-plane for NekRS-MOOSE, baseline THM-MOOSE, and the difference between the two (NekRS minus baseline THM).

Given the limitations of the 1-D area-averaged equations in THM, the temperature differences observed in Figure 6 and Figure 7 suggest that some differences between NekRS and THM could reasonably be ascribed to differences in spatial resolution and small differences in Nusselt number. A common multiscale use case of NekRS is the generation of closures for coarse-mesh thermal-fluid tools. In the next section, we will describe how Cardinal can be used to generate these heat transfer coefficients.

Heat Transfer Coefficients

In this section, we use Cardinal to correct the THM closures to better reflect the unit cell geometry and thermal-fluid conditions of interest. A correction is only computed at the single Reynolds and Prandtl numbers characterizing the unit cell ( and ). That is, a new coefficient is computed for a Dittus-Boelter type correlation with the same functional dependence on and ,

(4)

The heat transfer coefficient was defined in Eq. (3). We simply re-run the NekRS case and activate the [UserObjects] and [VectorPostprocessors] blocks (to save time, you can restart your NekRS case from one of the output files generated in the first part of this tutorial, making sure to comment out the temperature IC in the ranstube.udf because you'd now be loading temperature). Once you re-run the inputs, you will have CSV files that contain , , and in each axial layer that you can write a simple Python script to calculate and convert to a Nusselt number.

From the results of this additional simulation, Figure 8 shows the Nusselt number as a function of distance from the channel inlet computed by NekRS as a dashed line. A single average Nusselt number from the NekRS simulation is obtained as an axial average over the entire NekRS solution. Also shown for comparison is the Nusselt number predicted by the Dittus-Boelter correlation. NekRS predicts an average Nusselt number of 336.35, whereas the Dittus-Boelter correlation predicts a Nusselt number of 369.14; that is, the average NekRS Nusselt number is 8.9% smaller than the Dittus-Boelter correlation, which agrees very well with the qualitative assessment surrounding the discussion on Figure 6 that the leading contributor to a difference between the solid temperature with a NekRS fluid coupling vs. a THM fluid coupling could be ascribed to a small difference in convective heat transfer coefficient.

Figure 8: Nusselt number as a function of distance from the inlet for NekRS and the Dittus-Boelter correlation.

The NekRS simulation suggests a new coefficient for Eq. (4). The Dittus-Boelter correlation is updated in the THM input file by simply changing the 0.023 coefficient in the Nusselt number correlation to 0.021.

[Materials]
  [Nu_mat]
    type = ADParsedMaterial
    block = channel
    # Dittus-Boelter
    expression = '0.023 * pow(Re, 0.8) * pow(Pr, 0.4)'
    property_name = 'Nu'
    material_property_names = 'Re Pr'
  []
[]
(tutorials/gas_compact_cht/thm.i)

Then, the THM-MOOSE simulations repeated with the NekRS-informed closure. Figure 9 compares the solid temperature between NekRS-MOOSE with the corrected THM-MOOSE simulation. With the corrected Nusselt number correlation, the maximum solid temperature difference between NekRS-MOOSE and corrected THM-MOOSE is only 10.3 K. Due to the use of an axially-constant Nusselt number in THM, a perfect agreement between NekRS-MOOSE and corrected THM-MOOSE cannot be expected. The solid temperature difference will tend to be largest in regions where the local Nusselt number (dashed red line in Figure 8) deviates most from the global average (solid red line in Figure 8). This effect is observed near the outlet of the domain, where the local Nusselt number is less than half of the average Nusselt number, causing the highest deviation between NekRS-MOOSE and corrected THM-MOOSE. On average, the magnitude of the temperature difference is only 5.3 K.

Figure 9: Solid temperature predicted for NekRS-MOOSE, corrected THM-MOOSE, and the difference between the two (NekRS minus corrected THM).

Figure 10 shows the solid temperature predictions on the axial midplane for the NekRS-MOOSE simulation, the corrected THM-MOOSE simulation, and the difference between the two (NekRS minus corrected THM). Changing the Nusselt number model in THM only affects the wall temperature applied as a boundary condition in THM; because THM solves 1-D equations, the radial variation in the error will still be present (although of a smaller magnitude due to the temperature shift induced by a change in wall temperature).

Figure 10: Solid temperature predicted along the axial midplane for NekRS-MOOSE, corrected THM-MOOSE, and the difference between the two (NekRS minus corrected THM).

Cardinal's seamless capabilities for generating coarse-mesh closures using NekRS has been demonstrated. Now, just to present a few remaining results (all obtained with the corrected THMN-MOOSE simulation). Figure 11 shows the fluid temperature predicted by NekRS-MOOSE and the corrected THM-MOOSE. Also shown for context is the solid temperature in a slice through the domain, on a different color scale than the fluid temperature. The NekRS simulation resolves the fluid thermal boundary layer, whereas the THM model uses the Dittus-Boelter correlation to represent the temperature difference between the heated wall and the bulk; therefore, the fluid temperature shown in Figure 11 for THM is the area-averaged fluid temperature. The wall temperature predicted by NekRS follows a similar distribution as the heat flux, peaking slightly downstream of the midplane due to the combined effects of convective heat transfer and the imposed sinusoidal power distribution.

Figure 11: Fluid temperature predicted for NekRS-MOOSE and corrected THM-MOOSE.

Finally, Figure 12 compares the solutions for radially-averaged temperatures along the flow direction. As already discussed, the solid temperatures match very well between NekRS-MOOSE and THM-MOOSE and match the expected behavior for a sinusoidal power distribution. The fluid bulk temperature increases along the flow direction, with a faster rate of increase where the power density is highest.

Figure 12: Radially-averaged temperatures for the NekRS-MOOSE and corrected THM-MOOSE simulations.

References

  1. R. A. Berry, J. W. Peterson, H. Zhang, R. C. Martineau, H. Zhao, L. Zou, and D. Andrš. RELAP-7 Theory Manual. Technical Report INL/EXT-14-31366, Idaho National Laboratory, 2014.[BibTeX]
  2. F.P. Incropera and D.P. Dewitt. Fundamentals of Heat and Mass Transfer, 7th ed. John Wiley & Sons, 2011.[BibTeX]
  3. INL. High-Temperature Gas-Cooled Test Reactor Point Design. Technical Report INL/EXT-16-38296, Idaho National Laboratory, 2016. URL: https://tinyurl.com/v9a4hym5.[BibTeX]
  4. A.J. Novak, D. Andrs, P. Shriwise, J. Fang, H. Yuan, D. Shaver, E. Merzari, P.K. Romano, and R.C. Martineau. Coupled Monte Carlo and Thermal-Fluid Modeling of High Temperature Gas Reactors Using Cardinal. Annals of Nuclear Energy, 177:109310, 2022. doi:10.1016/j.anucene.2022.109310.[BibTeX]
  5. R.H.S. Winterton. Where did the Dittus and Boelter Equation Come From. International Journal of Heat and Mass Transfer, 41:809–810, 1998. doi:10.1016/S0017-9310(97)00177-4.[BibTeX]