Merge branch 'sfr_7pin' into 'master'

Sfr 7pin example

See merge request !9

sfr_7pin/duct.jou

+#!python
+
+# Model of a single wire pitch of the ABTR. All design specifications are from
+#
+# "Advanced Burner Test Reactor Preconceptual Design Report." ANL (2006)
+# https://tinyurl.com/ri5ps6jw
+
+# All units are SI (meters for length), corresponding to cold conditions
+# before the fuel-clad gap has closed.
+
+import math
+import os
+
+# append path to SFR mesh builder module
+import sys
+sys.path.append('/Users/anovak/projects/inputs/utils')
+import sfr_mesh as sfr_mesh
+
+cubit.cmd('reset')
+
+clad_OD = 8.0e-3
+P_D = 1.13
+pin_pitch = P_D * clad_OD
+n_pin_rings = 2
+length = 20.32e-2
+t_duct = 0.3e-2
+
+# the bundle flat-to-flat data is given for the 217-pin bundle, so we need to figure
+# out what the clearance is for the reduced 7-pin version
+full_flat_to_flat = 14.198e-2 - 2.0 * 0.3e-2
+full_bundle_flat_to_flat = sfr_mesh.min_pitch(sfr_mesh.rings(217), pin_pitch, clad_OD)
+pin_duct_clearance = full_flat_to_flat - full_bundle_flat_to_flat
+bundle_flat_to_flat = sfr_mesh.min_pitch(n_pin_rings, pin_pitch, clad_OD) + pin_duct_clearance
+
+fluid = sfr_mesh.HexagonalPrism(length, bundle_flat_to_flat)
+duct = sfr_mesh.HexagonalPrism(length, bundle_flat_to_flat + 2 * t_duct)
+
+# create the fluid region
+cubit.cmd(fluid.hexagonal_prism(1))
+
+# create the duct region
+cubit.cmd(duct.hexagonal_prism(2))
+
+cubit.cmd('subtract volume 1 from volume 2')
+
+# compress the surface IDs to begin at 1
+cubit.cmd('compress all')
+
+cubit.cmd(' move Volume 1 x 0 y 0 z ' + str(length / 2.0) + ' include_merged')
+
+# set the number of elements that will appear on each curve
+vertical_slices = 40
+cubit.cmd('curve 3 4 7 10 13 16 19 21 23 26 29 32 interval ' + str(vertical_slices))
+cubit.cmd('curve 3 4 7 10 13 16 19 21 23 26 29 32 scheme equal')
+
+# sides of the prism - horizontal curves; the azimuthal_slices is the number of elements
+# in the azimuthal direction. For symmetry purposes, this should be evenly divisible by 6
+wall_interval = 5
+
+cubit.cmd('curve 2 6 9 12 15 18 1 5 8 11 14 17 20 24 27 30 33 35 22 25 28 31 34 36 interval ' + str(math.ceil(wall_interval)))
+cubit.cmd('curve 2 6 9 12 15 18 1 5 8 11 14 17 20 24 27 30 33 35 22 25 28 31 34 36 scheme equal')
+
+cubit.cmd('mesh surface 1 2 3 4 5 6 7 8 9 10 11 12')
+
+# boundary layer on the hexagonal pincell surface
+cubit.cmd('create boundary_layer 1')
+cubit.cmd('modify boundary_layer 1 uniform height 4e-4 growth 1.2 layers 2')
+cubit.cmd('modify boundary_layer 1 add surface 1 volume 1 surface 2 volume 1 surface 3 volume 1 surface 4 volume 1 surface 5 volume 1 surface 6 volume 1')
+
+cubit.cmd('set boundary_layer intersection volume 1 curve 3 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 4 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 7 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 10 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 13 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 16 type side')
+
+cubit.cmd('create boundary_layer 2')
+cubit.cmd('modify boundary_layer 2 uniform height 4e-4 growth 1.2 layers 2')
+cubit.cmd('modify boundary_layer 2 add surface 7 volume 1 surface 8 volume 1 surface 9 volume 1 surface 10 volume 1 surface 11 volume 1 surface 12 volume 1')
+cubit.cmd('set boundary_layer intersection volume 1 curve 19 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 21 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 23 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 26 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 29 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 32 type side')
+
+cubit.cmd('mesh volume 1')
+
+# define volumes
+cubit.cmd('block 10 volume 1')
+cubit.cmd('block 10 element type HEX8')
+
+# define sidesets for applying BCs
+cubit.cmd('sideset 10 surface 7 8 9 10 11 12')
+cubit.cmd('sideset 10 name "duct_inner"')
+cubit.cmd('sideset 20 surface 1 2 3 4 5 6')
+cubit.cmd('sideset 20 name "duct_outer"')
+
+HOME = os.getenv('HOME')
+
+cubit.cmd('set large exodus file on')
+cubit.cmd('export Genesis "' + HOME + '/projects/inputs/cardinal/sfr_7pin/without_wire/duct.exo" dimension 3 overwrite')

sfr_7pin/mesh.jou

+#!python
+
+# Model of a single wire pitch of the ABTR. All design specifications are from
+#
+# "Advanced Burner Test Reactor Preconceptual Design Report." ANL (2006)
+# https://tinyurl.com/ri5ps6jw
+
+# All units are SI (meters for length), corresponding to cold conditions
+# before the fuel-clad gap has closed.
+
+import math
+import os
+
+# append path to SFR mesh builder module
+import sys
+sys.path.append('/Users/anovak/projects/inputs/utils')
+import sfr_mesh as sfr_mesh
+
+cubit.cmd('reset')
+
+clad_OD = 8.0e-3
+P_D = 1.13
+pin_pitch = P_D * clad_OD
+n_pin_rings = 2
+length = 20.32e-2
+
+# the bundle flat-to-flat data is given for the 217-pin bundle, so we need to figure
+# out what the clearance is for the reduced 7-pin version
+full_flat_to_flat = 14.198e-2 - 2.0 * 0.3e-2
+full_bundle_flat_to_flat = sfr_mesh.min_pitch(sfr_mesh.rings(217), pin_pitch, clad_OD)
+pin_duct_clearance = full_flat_to_flat - full_bundle_flat_to_flat
+bundle_flat_to_flat = sfr_mesh.min_pitch(n_pin_rings, pin_pitch, clad_OD) + pin_duct_clearance
+
+fluid = sfr_mesh.HexagonalPrism(length, bundle_flat_to_flat)
+
+# create the fluid region
+cubit.cmd(fluid.hexagonal_prism(1))
+
+# get the pin centers
+pin_centers = sfr_mesh.pin_centers(n_pin_rings, pin_pitch) 
+
+# create the pins
+n_pins = sfr_mesh.elements(n_pin_rings)
+for i in range(n_pins):
+  vol_id = i + 2
+  cubit.cmd('create cylinder height ' + str(length) + ' radius ' + str(clad_OD / 2.0))
+  cubit.cmd('move volume ' + str(vol_id) + ' location ' + str(pin_centers[i, 0]) + ' ' + str(pin_centers[i, 1]) + ' 0 include_merged')
+
+pin_volumes = ''
+for i in range(n_pins):
+  vol_id = i + 2
+  pin_volumes += str(vol_id) + ' '
+
+# get the fluid region only
+cubit.cmd('subtract volume ' + pin_volumes + ' from volume 1')
+
+# compress the surface IDs to begin at 1
+cubit.cmd('compress all')
+
+cubit.cmd(' move Volume 1 x 0 y 0 z ' + str(length / 2.0) + ' include_merged')
+
+# set the number of elements that will appear on each curve; to ensure a symmetric
+# transition from the cylinder surface to the hex walls, set the cylinder interval to
+# be six times the interval of each wall.
+
+# sides of the prism - vertical curves; the vertical_slices is the number of layers
+# desired in the axial direction
+vertical_slices = 40
+cubit.cmd('curve 3 4 7 10 13 16 interval ' + str(vertical_slices))
+cubit.cmd('curve 3 4 7 10 13 16 scheme equal')
+
+# sides of the prism - horizontal curves; the azimuthal_slices is the number of elements
+# in the azimuthal direction. For symmetry purposes, this should be evenly divisible by 6
+azimuthal_slices = 60
+arc_length_azimuthal = math.pi * clad_OD / azimuthal_slices
+
+wall_interval = sfr_mesh.side_length(bundle_flat_to_flat) / arc_length_azimuthal
+
+cubit.cmd('curve 2 6 9 12 15 18 interval ' + str(math.ceil(wall_interval)))
+cubit.cmd('curve 2 6 9 12 15 18 scheme equal')
+
+pin_curves = ''
+for i in range(n_pins):
+  curve_id = 2*i + 19
+  pin_curves += str(curve_id) + ' ' + str(curve_id + 1) + ' '
+  
+cubit.cmd('curve ' + pin_curves + ' interval ' + str(int(azimuthal_slices)))
+
+pin_surfaces = ''
+for i in range(n_pins):
+  surface_id = i + 7
+  pin_surfaces += str(surface_id) + ' '
+
+cubit.cmd('mesh surface 1 2 3 4 5 6')
+
+# boundary layer on the hexagonal pincell surface
+cubit.cmd('create boundary_layer 1')
+cubit.cmd('modify boundary_layer 1 uniform height 1e-4 growth 1.2 layers 3')
+cubit.cmd('modify boundary_layer 1 add surface 1 volume 1 surface 2 volume 1 surface 3 volume 1 surface 4 volume 1 surface 5 volume 1 surface 6 volume 1')
+
+cubit.cmd('set boundary_layer intersection volume 1 curve 3 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 4 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 7 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 10 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 13 type side')
+cubit.cmd('set boundary_layer intersection volume 1 curve 16 type side')
+
+# boundary layer on the pin
+cubit.cmd('create boundary_layer 2')
+cubit.cmd('modify boundary_layer 2 uniform height 1e-4 growth 1.2 layers 3')
+
+cmd = ''
+for i in range(n_pins):
+  surface_id = i + 7
+  cmd += 'surface ' + str(surface_id) + ' volume 1 '
+
+cubit.cmd('modify boundary_layer 2 add ' + cmd)
+
+# define volumes
+cubit.cmd('block 1 volume 1')
+
+cubit.cmd('mesh vol 1')
+cubit.cmd('block 1 element type HEX20')
+
+# define sidesets for applying BCs
+cubit.cmd('sideset 1 surface ' + pin_surfaces)
+cubit.cmd('sideset 2 surface 1 2 3 4 5 6')
+cubit.cmd('sideset 3 surface 15')
+cubit.cmd('sideset 4 surface 14')
+
+HOME = os.getenv('HOME')
+
+cubit.cmd('set large exodus file on')
+cubit.cmd('export Genesis "' + HOME + '/projects/inputs/cardinal/sfr_7pin/without_wire/sfr_7pin.exo" dimension 3 overwrite')
+
+print('Hydraulic diameter: ' + str(sfr_mesh.hydraulic_diameter(1, pin_pitch, clad_OD / 2.0)))
+print('Flow area: ' + str(sfr_mesh.flow_area(1, pin_pitch, clad_OD / 2.0)))
+
+for i in range(n_pins):
+ print('%+.8f %+.8f %+.8f' %(pin_centers[i, 0], pin_centers[i, 1], 0.0))

sfr_7pin/nek.i

+# This input file controls the nekRS run (in addition to the usual nekRS input files).
+#
+# -- sideset 1: pellet surface
+# -- sideset 2: duct inner surface
+# -- sideset 3: inlet
+# -- sideset 4: outlet
+
+# Note that these sideset IDs are different from what you see in solid.i! This is
+# because nekRS uses a completely separate mesh from MOOSE, and it so happens to have
+# different sideset IDs.
+
+[Mesh]
+  # First, build a mirror of the nekRS mesh, but only on the boundaries that will be
+  # communicating with MOOSE. Here, we want to send temperatures from nekRS and heat
+  # flux into nekRS on the pin surfaces and duct inner surface.
+  type = NekRSMesh
+  boundary = '1 2'
+[]
+
+[Problem]
+  type = NekRSProblem
+[]
+
+[Executioner]
+  type = Transient
+
+  [./TimeStepper]
+    type = NekTimeStepper
+  [../]
+[]
+
+[Outputs]
+  exodus = true
+  execute_on = 'final'
+[]
+
+[Postprocessors]
+  [flux_integral] # this receives the total heat flux for normalization
+    type = Receiver
+  []
+
+  # 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'
+  [pin_flux_in_nek]
+    type = NekHeatFluxIntegral
+    boundary = '1'
+  []
+  [duct_flux_in_nek]
+    type = NekHeatFluxIntegral
+    boundary = '2'
+  []
+
+  [max_nek_T]
+    type = NekVolumeExtremeValue
+    field = temperature
+    value_type = max
+  []
+  [min_nek_T]
+    type = NekVolumeExtremeValue
+    field = temperature
+    value_type = min
+  []
+[]

sfr_7pin/nek_master.i

+[Mesh]
+  # make a mesh of a single pincell
+  [circle]
+    type = AnnularMeshGenerator
+    nr = 20
+    nt = 20
+    rmin = 0
+    rmax = 0.0033
+    growth_r = -1.3
+  []
+  [cylinder]
+    type = FancyExtruderGenerator
+    input = circle
+    heights = '0.150'
+    num_layers = '20'
+    direction = '0 0 1'
+  []
+
+  # and then translate it to all the locations where there are pins
+  [combine]
+    type = CombinerGenerator
+    inputs = cylinder
+    positions = '
+      +0.00000000 +0.00000000 +0.00000000
+      +0.00414000 +0.00717069 +0.00000000
+      -0.00414000 +0.00717069 +0.00000000
+      -0.00828000 +0.00000000 +0.00000000
+      -0.00414000 -0.00717069 +0.00000000
+      +0.00414000 -0.00717069 +0.00000000
+      +0.00828000 +0.00000000 +0.00000000'
+  []
+
+  # make sure the mesh lines up with nekRS mesh, which is centered on (0, 0, 0)
+  [translate]
+    type = TransformGenerator
+    transform = translate
+    vector_value = '0 0 -0.075'
+    input = combine
+  []
+[]
+
+[Variables]
+  [temperature]
+    initial_condition = 500.0
+  []
+[]
+
+[Kernels]
+  [diffusion]
+    type = HeatConduction
+    variable = temperature
+    diffusion_coefficient = thermal_conductivity
+  []
+  [source]
+    type = CoupledForce
+    variable = temperature
+    v = source
+  []
+[]
+
+[BCs]
+  [interface]
+    type = MatchedValueBC
+    variable = temperature
+    v = nek_temp
+    boundary = '1'
+  []
+[]
+
+[Executioner]
+  type = Transient
+  num_steps = 1
+  dt = 0.05
+  nl_abs_tol = 1e-8
+[]
+
+[MultiApps]
+  [nek]
+    type = TransientMultiApp
+    app_type = NekApp
+    input_files = 'nek.i'
+    sub_cycling = true
+    execute_on = timestep_end
+  []
+[]
+
+[Transfers]
+  [temperature]
+    type = MultiAppNearestNodeTransfer
+    source_variable = temp
+    direction = from_multiapp
+    multi_app = nek
+    variable = nek_temp
+  []
+  [flux]
+    type = MultiAppNearestNodeTransfer
+    source_variable = flux
+    direction = to_multiapp
+    multi_app = nek
+    variable = avg_flux
+    source_boundary = '1'
+  []
+  [flux_integral]
+    type = MultiAppPostprocessorTransfer
+    to_postprocessor = flux_integral
+    direction = to_multiapp
+    from_postprocessor = flux_integral
+    multi_app = nek
+  []
+  [synchronization_to_nek]
+    type = MultiAppPostprocessorTransfer
+    direction = to_multiapp
+    to_postprocessor = synchronization_in
+    from_postprocessor = synchronization_to_nek
+    multi_app = nek
+  []
+[]
+
+[AuxVariables]
+  [flux]
+    family = MONOMIAL
+    order = CONSTANT
+  []
+  [nek_temp]
+    initial_condition = 628.15
+  []
+  [source]
+    initial_condition = 1e9
+  []
+[]
+
+[AuxKernels]
+  [flux]
+    type = NormalDiffusionFluxAux
+    variable = flux
+    coupled = temperature
+    diffusivity = thermal_conductivity
+    boundary = '1'
+  []
+[]
+
+[Materials]
+  [k]
+    type = GenericConstantMaterial
+    prop_names = 'thermal_conductivity'
+    prop_values = '1.5'
+  []
+[]
+
+[Postprocessors]
+active = 'flux_integral synchronization_to_nek'
+  [flux_integral]
+    type = SideIntegralVariablePostprocessor
+    variable = flux
+    boundary = '1'
+  []
+  [max_temp_bison]
+    type = NodalExtremeValue
+    variable = temperature
+    value_type = max
+  []
+  [min_temp_bison]
+    type = NodalExtremeValue
+    variable = temperature
+    value_type = min
+  []
+  [max_temp_interface]
+    type = NodalExtremeValue
+    variable = nek_temp
+    value_type = max
+  []
+  [min_temp_interface]
+    type = NodalExtremeValue
+    variable = nek_temp
+    value_type = min
+  []
+  [power]
+    type = ElementIntegralVariablePostprocessor
+    variable = source
+  []
+  [synchronization_to_nek]
+    type = Receiver
+    default = 1
+  []
+[]
+
+[Outputs]
+  exodus = true
+  print_linear_residuals = false
+[]

sfr_7pin/sfr_7pin.oudf

+// This is almost a normal .oudf file - we just need to be sure that for the
+// boundaries that will receive heat flux from MOOSE, to set it with
+// bc->wrk[bc->idM] (we write into this array from MOOSE).
+
+void velocityDirichletConditions(bcData * bc)
+{
+  switch (bc->id)
+  {
+    case 3:
+      bc->u = 0.0; // x-velocity
+      bc->v = 0.0; // y-velocity
+      bc->w = Vz;  // z-velocity
+      break;
+    default:
+       throw std::runtime_error("Invalid Dirichlet velocity BC sideset! You entered: " +
+         std::to_string(bc->id) + "; valid options: 3.");
+  }
+}
+
+void scalarDirichletConditions(bcData * bc)
+{
+  if (bc->scalarId == 0) // temperature
+  {
+    switch (bc->id)
+    {
+      case 3:
+        bc->s = inlet_T;
+        break;
+      default:
+         throw std::runtime_error("Invalid Dirichlet temperature BC sideset! You entered: " +
+           std::to_string(bc->id) + "; valid options: 3.");
+    }
+  }
+  else
+  {
+    throw std::runtime_error("Invalid scalar ID for BC application! You entered: " +
+      std::to_string(bc->scalarId) + "; valid options: 0.");
+  }
+}
+
+void scalarNeumannConditions(bcData * bc)
+{
+  if (bc->scalarId == 0) // temperature
+  {
+    switch (bc->id)
+    {
+      case 1: // boundary coupled to MOOSE
+        bc->flux = bc->wrk[bc->idM];
+        break;
+      case 2: // boundary coupled to MOOSE
+        bc->flux = bc->wrk[bc->idM];
+        break;
+      default:
+        break;
+         throw std::runtime_error("Invalid Neumann temperature BC sideset! You entered: " +
+           std::to_string(bc->id) + "; valid options: 1, 2.");
+    }
+  }
+  else
+  {
+    throw std::runtime_error("Invalid scalar ID for BC application! You entered: " +
+      std::to_string(bc->scalarId) + "; valid options: 0.");
+  }
+}

sfr_7pin/sfr_7pin.par

+[GENERAL]
+  stopAt = numSteps
+  numSteps = 1000
+  dt = 5.0e-4
+  variableDT = no
+  targetCFL = 0.5
+  timeStepper = tombo2
+  extrapolation = standard
+  subCyclingSteps = 1
+  writeControl = timeStep
+  writeInterval = 25
+  filtering = none
+  dealiasing = yes
+  polynomialOrder = 3
+  optLevel = 2
+
+[PROBLEMTYPE]
+  equation = incompNS
+  axiSymmetry = no
+  swirl = no
+  cyclicBoundaries = no
+  variableProperties = yes
+  stressFormulation = yes
+
+[MESH]
+  motion = none
+
+[VELOCITY]
+  viscosity = 2.37e-4
+  density = 834.5
+  boundaryTypeMap = W, W, v, O
+  residualTol = 1.0e-6
+  residualProj = no
+
+[PRESSURE]
+  residualTol = 1.0e-6
+  residualProj = yes
+  preconditioner = semg_amg
+
+[TEMPERATURE]
+  conductivity = 64.21
+  rhoCp = 1024766.0
+
+  # the first two boundaries are where we will exchange data with MOOSE - we receive
+  # a heat flux from MOOSE, so we need to set these to flux BCs
+  boundaryTypeMap = f, f, t, O
+  residualTol = 1.0e-5
+  absoluteTol = 1.0e-5

sfr_7pin/sfr_7pin.udf

+#include "udf.hpp"
+#include "bcMap.hpp"
+
+// This is just a usual UDF file - nothing special for coupling here
+
+float mass_flowrate()
+{
+  // TODO: add turbulence and then increase this to a realistic value
+  return 0.1;
+}
+
+float height()
+{
+  return 20.32e-2;
+}
+
+float inlet_temp()
+{
+  return 355 + 273.15; // page 25
+}
+
+float inlet_velocity(nrs_t * nrs)
+{
+  setupAide & options = nrs->options;
+
+  // flow area (comes from Cubit area measurement)
+  float flow_area = 0.000245;
+
+  // get density from input file
+  float rho;
+  options.getArgs("DENSITY", rho);
+
+  return mass_flowrate() / flow_area / rho;
+}
+
+void UDF_LoadKernels(nrs_t *nrs)
+{
+  setupAide & options = nrs->options;
+  occa::properties & kernelInfo = *nrs->kernelInfo;
+
+  // inlet axial velocity (assuming zero radial and azimuthal components)
+  kernelInfo["defines/Vz"] = inlet_velocity(nrs);
+
+  // inlet temperature
+  kernelInfo["defines/inlet_T"] = inlet_temp();
+}
+
+void UDF_Setup0(MPI_Comm comm, setupAide &options)
+{
+}
+
+void UDF_Setup(nrs_t *nrs)
+{
+  // set initial conditions for the velocity, temperature, and pressure
+  mesh_t * mesh = nrs->mesh;
+
+  float inlet_vel = inlet_velocity(nrs);
+
+  // loop over all the GLL points and assign directly to the solution arrays by
+  // indexing according to the field offset necessary to hold the data for each
+  // solution component
+  int n_gll_points = mesh->Np * mesh->Nelements;
+  for (int n = 0; n < n_gll_points; ++n)
+  {
+    nrs->U[n + 0 * nrs->fieldOffset] = 0.0;
+    nrs->U[n + 1 * nrs->fieldOffset] = 0.0;
+    nrs->U[n + 2 * nrs->fieldOffset] = inlet_vel;
+
+    nrs->P[n] = 0.0;
+
+    nrs->cds->S[n + 0 * nrs->cds->fieldOffset] = inlet_temp();
+  }
+}