Commit 093f9550 authored by April Novak's avatar April Novak

Add example of a 7-pin bundle with explanations of the input files.

parent b4b0bd5b
#!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')
#!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))
# 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
[]
[]
[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
[]
// 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.");
}
}
[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
#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();
}
}
This diff is collapsed.
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment