Conjugate Heat Transfer Simulation of a 67-Pebble Core

Contact: John Acierno (Penn State), Dillon Shaver (dshaver.at.anl.gov), and April Novak (anovak.at.anl.gov)

Model link: 67-Pebble Core Model

In this tutorial we are going to set up and simulate a simple Conjugate Heat Transfer (CHT) case using a helium (Pr=0.71) cooled 67-pebble bed. This tutorial was developed from an example case provided with NekRS and couples to MOOSE's CHT module using Cardinal as a wrapper. More information about the NEAMS tool Cardinal can be found on github, or on the Cardinal website. In each time step MOOSE will solve the energy equation in the the solid subdomain and pass the solution to NekRS, which will in-turn solve both the Navier-Stokes and energy equations in the fluid subdomain. NekRS will then pass its temperature solution back to MOOSE in the next time step. This transfer of information occurs at the boundary between solid and fluid subdomains, which are the pebble surfaces in this case.

Cardinal and MOOSE-Wrapped Apps

The analysis shown here is performed with Cardinal, a MOOSE application that "wraps" nekRS, a Computational Fluid Dynamics (CFD) code, and OpenMC, a Monte Carlo particle transport code. Cardinal is an open-source application, available on GitHub. As this example focuses on heat transfer modeling, there will be no further discussion of the OpenMC wrapping in Cardinal. "Wrapping" means that, for all intents and purposes, nekRS simulations can be run within the MOOSE framework and interacted with as if the physics and numerical solution performed by nekRS were a native MOOSE application. Cardinal contains source code that facilitates data transfers in and out of nekRS and runs nekRS within a MOOSE-controlled simulation. At a high level, Cardinal's wrapping of nekRS consists of:

  • Constructing a "mirror" of the nekRS mesh through which data transfers occur with MOOSE. For conjugate heat transfer applications such as those shown here, a MooseMesh is created by copying the nekRS surface mesh into a format that all native MOOSE applications can understand.

  • Adding MooseVariables to represent the nekRS solution. In other words, if nekRS stores the temperature internally as an std::vector<double>, with each entry corresponding to a nekRS node, then a MooseVariable is created that represents the same data, but that can be accessed in relation to the MooseMesh mirror.

  • Writing multiphysics feedback fields in/out of nekRS's internal solution and boundary condition arrays. So, if nekRS represents a heat flux boundary condition internally as an std::vector<double>, this involves reading from a MooseVariable representing heat flux (which can be transferred with any of MOOSE's transfers to the nekRS wrapping) and writing into nekRS's internal vectors.

Accomplishing the above three tasks requires an intimate knowledge of how nekRS stores its solution fields and mesh. But once the wrapping is constructed, nekRS can then communicate with any other MOOSE application via the MultiApp and Transfer systems in MOOSE, enabling complex multiscale thermal-hydraulic analysis and multiphysics feedback. Note that the user does not need to construct this mapping themselves - it all happens automatically. The same wrapping can be used for conjugate heat transfer analysis with any MOOSE application that can compute a heat flux; that is, because a MOOSE-wrapped version of nekRS interacts with the MOOSE framework in a similar manner as natively-developed MOOSE applications, the agnostic formulations of the MultiApps and Transfers can be used to equally extract heat flux from Pronghorn, BISON, the MOOSE heat conduction module, and so on.

schooltip:Why Use Cardinal?

Cardinal allows nekRS and OpenMC to couple seamlessly with the MOOSE framework, enabling in-memory coupling, distributed parallel meshes for very large-scale applications, and straightforward multiphysics problem setup. Cardinal has capabilities for conjugate heat transfer coupling of nekRS and MOOSE, concurrent conjugate heat transfer and volumetric heat source coupling of nekRS and MOOSE, and volumetric density, temperature, and heat source coupling of OpenMC to MOOSE. Together, the OpenMC and nekRS wrappings augment MOOSE by expanding the framework to include high-resolution spectral element CFD and Monte Carlo particle transport. For conjugate heat transfer applications in particular, the libMesh-based interpolations of fields between meshes enables fluid-solid heat transfer simulations on meshes that are not necessarily continuous on phase interfaces, allowing mesh resolution to be specified based on the underlying physics, rather than rigid continuity restrictions in single-application heat transfer codes.

NekRS

nekRS solves the Navier-Stokes and energy conservation equations. They can be solved in either low-Mach or incompressible forms. This model uses the incompressible form,

∂ui∂xi=0ρf(∂ui∂t+uj∂ui∂xj)=−∂P∂xi+∂τij∂xj+ρffiρfCp,f(∂Tf∂t+ui∂Tf∂xi)=∂∂xi(kf∂Tf∂xi)+q˙f \begin{aligned} \frac{\partial u_i}{\partial x_i} & =0 \\ \rho_f\left(\frac{\partial u_i}{\partial t}+ u_j\frac{\partial u_i}{\partial x_j}\right) & = -\frac{\partial P}{\partial x_i}+\frac{\partial \tau_{ij}}{\partial x_j} +\rho_f f_i \\ \rho_f C_{p,f}\left(\frac{\partial T_f}{\partial t}+ u_i\frac{\partial T_f}{\partial x_i}\right) & = \frac{\partial }{\partial x_i}\left( k_f\frac{\partial T_f}{\partial x_i}\right)+\dot{q}_f \end{aligned}(1)

where uiu_i is the fluid velocity, ρf\rho_f is the fluid density, PP is the pressure, τij\tau_{ij} is the viscous stress tensor, fif_i is a body force vector, q˙f\dot{q}_f is a volumetric heat source in the fluid domain, Cp,fC_{p,f} is the fluid heat capacity, TfT_f is the fluid temperature, and kfk_f is the fluid thermal conductivity.

commentnote

The form of the viscous stress tensor depends on the problem formulation and the selected turbulence model.

nekRS can be solved in either dimensional or non-dimensional form; in the present case, nekRS is solved in non-dimensional form. That is, characteristic scales for velocity, temperature, length, and time are defined and substituted into the governing equations so that all solution fields (velocity, pressure, temperature) are of order unity. A full derivation of the non-dimensional governing equations in nekRS is available here. To summarize, the non-dimensional quantities are defined as

ui†≡uiU0xi†≡xiL0t†≡tU0L0P†≡PρfU02fi†≡fiL0U02T†≡T−T0ΔT0q˙†≡q˙fL0ρfU0Cp,fΔT0\begin{aligned} u_i^\dagger & \equiv \frac{u_i}{U_{0}} \\ x_i^\dagger & \equiv \frac{x_i}{L_{0}} \\ t^\dagger & \equiv t \frac{U_{0}}{L_{0}} \\ P^\dagger & \equiv \frac{P}{\rho_f U_{0}^2} \\ f_i^\dagger & \equiv f_i \frac{L_{0}}{U_{0}^2} \\ T^\dagger & \equiv \frac{T-T_{0}}{\Delta T_{0}} \\ \dot q^\dagger & \equiv \dot q_f \frac{L_{0}}{\rho_f U_{0} C_{p,f}\Delta T_{0}} \end{aligned}

where a †\dagger superscript indicates a non-dimensional quantity and a "00" subscript indicates a characteristic value. Inserting these definitions into Eq. (1) gives

∂ui†∂xi†=0∂ui†∂t†+uj†∂ui†∂xj†=−∂P†∂xi†+1Re∂τij†∂xj†+fi†∂T†∂t†+ui†∂T†∂xi†=1Pe∂2T†∂xi†2+q˙† \begin{aligned} \frac{\partial u_i^\dagger}{\partial x_i^\dagger} & = 0 \\ \frac{\partial u_i^\dagger}{\partial t^\dagger}+u_j^\dagger\frac{\partial u_i^\dagger}{\partial x_j^\dagger} & =-\frac{\partial P^\dagger}{\partial x_i^\dagger}+\frac{1}{Re}\frac{\partial \tau_{ij}^\dagger}{\partial x_j^\dagger}+f_i^\dagger \\ \frac{\partial T^\dagger}{\partial t^\dagger}+u_i^\dagger\frac{\partial T^\dagger}{\partial x_i^\dagger} & =\frac{1}{Pe}\frac{\partial^2 T^\dagger}{\partial x_i^{\dagger 2}}+\dot{q}^\dagger \end{aligned}(2)

Complete definitions of all terms are available here. New terms in these non-dimensional equations are ReRe and PePe, the Reynolds and Peclet numbers, respectively:

Re≡ρfU0L0μfPe≡ρfU0L0Cp,fkf=RePr\begin{aligned} Re & \equiv\frac{\rho_f U_{0}L_{0}}{\mu_f} \\ Pe & \equiv\frac{\rho_f U_{0} L_{0} C_{p,f}}{k_f}=RePr \end{aligned}

nekRS solves for ui†u_i^\dagger, P†P^\dagger, and T†T^\dagger. Cardinal will handle conversions from a non-dimensional solution to a dimensional MOOSE heat conduction application, but awareness of this non-dimensional formulation will be important for describing the nekRS input files.

This case uses an LES turbulence model in nekRS. For LES, the viscous stress tensor is given by

τij=μf(∂ui∂xj+∂uj∂xi)\tau_{ij} = \mu_f\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)

and a filtering operation is implemented as a source term to the momentum equation. More details about Nek filters can be found in Nek5000 Online Tutorial.

fi=−χ(ui−H∗ui)+f^if_i = - \chi \left(u_i - H * u_i\right) + \hat f_i

where χ\chi is the filter strength and (H ∗)(H\ *) is an element-wise convolution operator that acts on the specified number of modes (usually two). These are integrated into the dimensionless momentum equation to yield

∂ui†∂t†+uj†∂ui†∂xj†=−∂P†∂xi†+1Re∂2ui†∂xj†2−χ(ui†−H∗ui†)+f^i† \frac{\partial u_i^\dagger}{\partial t^\dagger}+u_j^\dagger\frac{\partial u_i^\dagger}{\partial x_j^\dagger} =-\frac{\partial P^\dagger}{\partial x_i^\dagger} +\frac{1}{Re}\frac{\partial^2 u_i^\dagger}{\partial x_j^{\dagger 2}} -\chi\left( u_i^\dagger - H * u_i^\dagger\right) +\hat f_i^\dagger(3)

where f^i†\hat f_i^\dagger may include other momentum source terms.

Heat Conduction Model

The MOOSE heat conduction module is used to solve for steady-state heat conduction,

−∇⋅(ks∇Ts)=q˙s -\nabla\cdot\left(k_s\nabla T_s\right)=\dot{q}_s(4)

where ksk_s is the solid thermal conductivity, TsT_s is the solid temperature, and q˙s\dot{q}_s is a volumetric heat source in the solid.

Geometry and computational model

The geometry for this case consists of 67 pebbles in a randomly packed cylindrical bed. The pebbles have a diameter of 1 (dimensionless) and the cylinder has a diameter of 4.4 (dimensionless). The pebble bed begins around 2.5 pebble diameters upstream from the cylinder inlet and has a total height of about 5 pebble diameters. A cross section of the fluid domain is shown in Figure 1 where the fluid is the gray region and the pebbles are the white regions.

Figure 1: A cross section of the fluid domain for the 67 pebble geometry.

NekRS Setup

For this case, the mesh file is provided through github Large File Storage (LFS) system. It can be downloaded in the following way

$ cd virtual_test_bed/htgr/pb67_cardinal
$ git lfs fetch pb67.re2

The mesh includes over 110,000 hexahedral elements and has been generated using a Voronoi cell approach Lan et al. (2021).

oudf file

We will begin by setting up the .oudf file. This file contains the setup for the required boundary conditions.

// Boundary conditions
void velocityDirichletConditions(bcData *bc)
{
  bc->u = 0.0;
  bc->v = 0.0;
  bc->w = 1.0;
}

// Stabilized outflow (Dong et al)
void pressureDirichletConditions(bcData *bc)
{
  const dfloat iU0delta = 20.0;
  const dfloat un = bc->u*bc->nx + bc->v*bc->ny + bc->w*bc->nz;
  const dfloat s0 = 0.5 * (1.0 - tanh(un*iU0delta));
  bc->p = -0.5 * (bc->u*bc->u + bc->v*bc->v + bc->w*bc->w) * s0;
}

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

void scalarNeumannConditions(bcData *bc)
{
  bc->flux = bc->usrwrk[bc->idM]; //get flux BC from MOOSE
}

@kernel void cliptOKL(const dlong Nelements,
                 @restrict dfloat * TEMP)
{
 for(dlong e=0;e<Nelements;++e;@outer(0)){
   for(int n=0;n<p_Np;++n;@inner(0)){
     const int id = e*p_Np + n;
     if(TEMP[id]>100.0)
          {
          TEMP[id] = 100.0;
          }
     if(TEMP[id]<0.0)
          {
          TEMP[id] = 0.0;
          }
  }
 }
}
(htgr/pb67_cardinal/pb67.oudf)

Note that we have bc->flux = bc->usrwrk[bc->idM] in scalarNeumannConditions block. This line allows NekRS to use the heat flux provided by the MOOSE CHT module at the surface of the pebbles instead of using a constant heat flux as previously defined. There is also a kernel function defined in the .oudf file, cliptOKL. The cliptOKL kernel is used to limit extreme temperatures in the simulation which can occur in underresolved parts of the mesh. If the temperature is greater than 100 or less than 0, this kernel will set the temperature to 100 or 0 respectively.

udf file

The .udf file utilizes the kernal functions we defined in .oudf, and provides the detailed configurations of the related nekRS simulations.

//
// nekRS User Defined File
//

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

static occa::kernel cliptKernel; // clipping

void clipt(nrs_t *nrs)
{
  mesh_t *mesh = nrs->meshV;
  cds_t* cds = nrs->cds;
  cliptKernel(mesh->Nelements, cds->o_S);
}

/* UDF Functions */

void UDF_LoadKernels(occa::properties & kernelInfo)
{
 // avg::buildKernel(nrs);
  cliptKernel = oudfBuildKernel(kernelInfo, "cliptOKL");
}

void UDF_Setup(nrs_t *nrs)
{
}

void UDF_ExecuteStep(nrs_t *nrs, dfloat time, int tstep)
{
  clipt(nrs);
  if (nrs->isOutputStep) {
    nek::ocopyToNek(time, tstep);
    nek::userchk();
  }
}
(htgr/pb67_cardinal/pb67.udf)

par file

The .par file lists the necessary input parameters to nekRS simulations.

[GENERAL]
polynomialOrder = 7
#startFrom = "restart.fld"
stopAt = endTime
endTime = 500

dt = 1.0e-05
timeStepper = tombo2
subCyclingSteps = 2

writeControl = steps
writeInterval = 500

filtering = hpfrt
filterWeight = 0.2/${dt}
filterModes = 2

[PRESSURE]
residualTol = 1e-04
preconditioner = multigrid
smootherType = chebyshev+jac
initialGuess = projectionAconj + nVector = 30

[VELOCITY]
boundaryTypeMap = inlet, outlet, wall, wall
density = 1.0
viscosity = -1460.0
residualTol = 1e-06

[TEMPERATURE]
boundaryTypeMap = t,I,I,f
residualTol = 1e-06
rhoCp = 1.0
conductivity = -1036.0
(htgr/pb67_cardinal/pb67.par)

If you want to run this example on GPUs, you simply need to pass --nekrs-backend CUDA --nekrs-device-id 0 on the command line later when we run the cardinal-opt executable. Otherwise, this example will be run on a CPU backend. The specific filtering settings allow us to smooth out high frequencies. If you are planning on doing a DNS you should not include any filtering. Next, the parameters in [PRESSURE] block are selected to have a short time-to-solution. We are able to achieve this by adding a preconditioner and a smoother. In NekRS version 21.1 the combination of semg preconditioner with chebyshev+asm smoothers works well, but in later versions 22.0+, using the amgx preconditioner with chebyshev+jac smoothers is much faster.

usr file

The .usr is a legacy from Nek5000 code, and we still use to specify the boundary conditons in Nek CFD calculations. Here we prescribe every element a boundaryID for both velocity and temperature. We are able to check the type of element using cbc(ifc,iel,1).eq.'TOP' and assign boundary types based on that. In this case we have 5 types of elements: inlet, outlet, outer wall, pebble wall, and chamfer wall. Notice, for the chamfer we assign the wall condition for velocity, but the E condition for temperature. This results in flow moving around the chamfered area without an additional heat flux. The chamfer is necessary to prevent a mesh singularity between touching pebbles.

subroutine usrdat2()  ! This routine to modify mesh coordinates
      include 'SIZE'
      include 'TOTAL'

      integer e,f

      do iel=1,nelt
      do ifc=1,2*ndim
      if (cbc(ifc,iel,1).eq.'TOP') then  ! top surface
          cbc(ifc,iel,1) = 'O  '
          cbc(ifc,iel,2) = 'I  '
          bc(5,ifc,iel,1) = 1
       else if  (cbc(ifc,iel,1).eq.'BOT') then  ! bot surface
          cbc(ifc,iel,1) = 'v  '
          cbc(ifc,iel,2) = 't  '
          bc(5,ifc,iel,1) = 2
       else if  (cbc(ifc,iel,1).eq.'SW ') then    ! side wall of cylinder
          cbc(ifc,iel,1) = 'W  '
          cbc(ifc,iel,2) = 'I  '
          bc(5,ifc,iel,1) = 3
       else if  (cbc(ifc,iel,1).eq.'PW ') then  ! pebble surface
          cbc(ifc,iel,1) = 'WH  '
          cbc(ifc,iel,2) = 'fH  '
          bc(5,ifc,iel,1) = 4
       else if  (cbc(ifc,iel,1).eq.'C  ') then  ! chamfer surface
          cbc(ifc,iel,1) = 'W  '
          cbc(ifc,iel,2) = 'E  '
          bc(5,ifc,iel,1) = 5
       else
          cbc(ifc,iel,1) = 'E  '
          bc(5,ifc,iel,1) = 0
      endif
      enddo
      enddo

      do iel=1,nelt
      do ifc=1,2*ndim
        boundaryID(ifc,iel)=0
        if(cbc(ifc,iel,1) .eq. 'v ') boundaryID(ifc,iel)=1
        if(cbc(ifc,iel,1) .eq. 'O ') boundaryID(ifc,iel)=2
        if(cbc(ifc,iel,1) .eq. 'W ') boundaryID(ifc,iel)=3
        if(cbc(ifc,iel,1) .eq. 'WH ') boundaryID(ifc,iel)=4
      enddo
      enddo

      do iel=1,nelt
      do ifc=1,2*ndim
         boundaryIDt(ifc,iel) = 0
         if (cbc(ifc,iel,2) .eq. 't  ') boundaryIDt(ifc,iel) = 1
         if (cbc(ifc,iel,2) .eq. 'I  ') boundaryIDt(ifc,iel) = 2
         if (cbc(ifc,iel,2) .eq. 'I  ') boundaryIDt(ifc,iel) = 3
         if (cbc(ifc,iel,2) .eq. 'fH  ') boundaryIDt(ifc,iel) = 4
      enddo
      enddo

      return
      end
c-----------------------------------------------------------------------
(htgr/pb67_cardinal/pb67.usr)

MOOSE setup

As previously stated, we strongly recommend completing the CHT tutorials on Cardinal, which fully describe every aspect of the MOOSE setup. Here, we will only cover case specific necessities.

In this case we need to create two files: nek.i and moose.i pertaining to the NekRS and MOOSE parameters respectively.

Starting with the simpler nek.i which contains the following

[Mesh]
  type = NekRSMesh
  boundary = 4
[]

[Problem]
  type = NekRSProblem
  casename = 'pb67'
[]

[Executioner]
  type = Transient
  [TimeStepper]
    type = NekTimeStepper
  []
[]

[Outputs]
  exodus = true
  interval=1000
[]

[Postprocessors]
  # This is the heat flux in the nekRS solution, i.e. it is not an integral
  # of nrs->usrwrk, instead this is directly an integral of k*grad(T)*hat(n).
  # So this should closely match 'flux_integral'
  [flux_in_nek]
    type = NekHeatFluxIntegral
    boundary = '4'
  []

  [max_nek_T]
    type = NekVolumeExtremeValue
    field = temperature
    value_type = max
  []
  [min_nek_T]
    type = NekVolumeExtremeValue
    field = temperature
    value_type = min
  []
  [average_nek_pebble_T]
    type = NekSideAverage
    boundary = '4'
    field = temperature
  []
[]
(htgr/pb67_cardinal/nek.i)

In the [Mesh] block, MOOSE will create a copy of the Nek mesh at the given boundaries. In this case it is at the pebble surfaces. In the [Problem] block we define the name of corresponding NekRS files. Next, we want MOOSE to use the same time steps as Nek so we prescribe that in the [Executioner] block. In the [Output] block we tell MOOSE to output an exodus file of the shallow Nek copy every 1,000 time steps. This can be a helpful check at the beginning of the simulation to make sure you are using the correct boundaries. Finally, in the [Postprocessor] block we define what values we want MOOSE to calculate for us. Here we want the integral flux, min, max, and average temperature at the pebble surface.

commentnote

Note that in this file, we are copying the non-dimensional nekRS solution into Cardinal. In this case, the coupled MOOSE heat conduction will also be represented in non-dimensional form. If we wanted to instead transform the non-dimensional nekRS solution into an equivalent dimensional form (which would allow us to represent the solid heat conduction with dimensional units), we'd simply set nondimensional = true for NekRSProblem and then set the characteristic scales, U_ref, T_ref, dT_ref, L_ref, rho_0, and Cp_0. For more information, see the NekRSProblem documentation.

To setup the MOOSE mesh, a text file containing the center points of each pebble is necessary. For the bed used in this case, the file can be obtained here.

[Mesh]
  [file]
    type = FileMeshGenerator
    file = sphere.e
  []
  [cmbn]
    type = CombinerGenerator
    inputs = file
    positions_file = 'positions.txt'
  []
  [scale]
    type = TransformGenerator
    input = cmbn
    transform = SCALE
    vector_value = '0.99 0.99 0.99'
  []
[]

[Kernels]
  [hc]
    type = HeatConduction
    variable = temp
  []
  [heat]
    type = BodyForce
    value = 0.01
    variable = temp
  []
[]

[BCs]
  [match_nek]
    type = MatchedValueBC
    variable = temp
    boundary = '1'
    v = 'nek_temp'
  []
[]

[Materials]
  [hc]
    type = GenericConstantMaterial
    prop_values = '0.005' # Should be 45 times higher than the fluid conductivity (Steel/P-cymene ratio)
    prop_names = 'thermal_conductivity'
  []
[]

[Executioner]
  type = Transient
  #  petsc_options_iname = '-pc_type -pc_hypre_type'
  num_steps = 16000
  #  petsc_options_value = 'hypre boomeramg'
  dt = 2e-4
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-10
[]

[Variables]
  [temp]
    initial_condition = 4.3
  []
[]

[Outputs]
  exodus = true
  time_step_interval = 500
  [csv]
    type = CSV
    time_step_interval = 1
  []
[]

[MultiApps]
  [nek]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'nek.i'
  []
[]

[Transfers]
  [nek_temp]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = temp
    from_multi_app = nek
    variable = nek_temp

    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [avg_flux]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = avg_flux
    to_multi_app = nek
    variable = avg_flux

    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [flux_integral_to_nek]
    type = MultiAppPostprocessorTransfer
    to_postprocessor = flux_integral
    from_postprocessor = flux_integral
    to_multi_app = nek
  []
[]

[AuxVariables]
  [nek_temp]
  []
  [avg_flux]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [avg_flux]
    type = DiffusionFluxAux
    diffusion_variable = temp
    component = normal
    diffusivity = thermal_conductivity
    variable = avg_flux
    boundary = '1'
  []
[]

[Postprocessors]
  [flux_integral]
    type = SideDiffusiveFluxIntegral
    diffusivity = thermal_conductivity
    variable = 'temp'
    boundary = '1'
  []
  [average_flux]
    type = SideDiffusiveFluxAverage
    diffusivity = thermal_conductivity
    variable = 'temp'
    boundary = '1'
  []
  [max_pebble_T]
    type = NodalExtremeValue
    variable = temp
    value_type = max
  []
  [min_pebble_T]
    type = NodalExtremeValue
    variable = temp
    value_type = min
  []
  [average_pebble_T]
    type = SideAverageValue
    variable = 'temp'
    boundary = '1'
  []
[]
(htgr/pb67_cardinal/moose.i)

Notice in the [Mesh] block we are using a file to generate the solid mesh. Make sure to include the sphere.e file in your working directory. We use this single sphere of radius 1 and duplicate it using CombinerGenerator and all of the coordinates of the 67 pebbles given in positions.txt. From here we scale down the pebbles to ensure that they are not touching in the solid mesh. We then call the heat conduction module in the [Kernels] block ensuring MOOSE will solve the temperature in the solid mesh. Then we tell MOOSE to match the boundary conditions use in Nek in the [BCs] block. In the [Transfers] block we define what and how variables get transfered. We define a multi app and transfer the temperature from Nek to MOOSE and return the average flux and flux integral back to Nek. Built in post-processors are used to obtain these values. Finally, we define an aux kernel to model the avg flux. Notice in nek.i the boundary from nek is 4 and in moose.i the boundary of interest is 1. Be sure not to mix up the boundaries in the fluid and solid mesh.

Results

Below are 3D renderings done in VisIt of the velocity (left) and temperature (right) fields.

Figure 2: 3D volumetric renderings done in VisIt of the velocity (left) and the temperature (right) fields.

References

  1. Yu-Hsiang Lan, Paul Fischer, Elia Merzari, and Misun Min. All-hex meshing strategies for densely packed spheres. arXiv, 2021. URL: https://arxiv.org/abs/2106.00196.[BibTeX]