LFR assembly - a problem with coupling MOOSE SubChannel and OpenMC

[albottar]

Good evening. There's a problem with coupling, or its just a bad design, I dont know. Im trying to model Heat Transfer - SCM - OpenMC system for LFR assembly (via tutorial (https://cardinal.cels.anl.gov/tutorials/subchannel.html)), but my dT about 10 K instead of 110 K with 5.5e6 Wth for assembly. Can anyone help me? I'm quite new to MOOSE, maybe there's a baby mistake somewhere.

My goal is to get a good temperature distribution for this model for future research.

`!include common.i

[Mesh]
  [file]
    type = FileMeshGenerator
    file = mesh_in.e
  []
  [delete_coolant]
    type = BlockDeletionGenerator
    input = file
    block = 'lead'
  []
[]

[Variables]
  [T]
  []
[]

[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = T
  []
  [heat_source_fuel]
    type = CoupledForce
    variable = T
    v = heat_source
    block = 'fuel'
  []
[]

[Materials]
   [k_helium]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity'
    prop_values = '1.5'
    block = 'helium'
  []
   [k_fuel]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity'
    prop_values = '26'
    block = 'fuel'
  []
   [k_clad]
    type = GenericConstantMaterial
    prop_names = 'thermal_conductivity'
    prop_values = '20'
    block = 'clad'
  []
   [fuel_density]
    type = GenericConstantMaterial
    prop_names = 'density'
    prop_values = '14600'
    block = 'fuel'
  []
   [helium_density]
   type = GenericConstantMaterial
    prop_names = 'density'
    prop_values = '1.7e-3'
    block = 'helium'
  []   
   [clad_density]
   type = GenericConstantMaterial
    prop_names = 'density'
    prop_values = '6000'
    block = 'clad'
  []
   [fuel_cp]
    type = GenericConstantMaterial
    prop_names = 'specific_heat'
    prop_values = '234'
    block = 'fuel'
  []
   [helium_cp]
   type = GenericConstantMaterial
    prop_names = 'specific_heat'
    prop_values = '3300'
    block = 'helium'
  []   
   [clad_cp]
   type = GenericConstantMaterial
    prop_names = 'specific_heat'
    prop_values = '291'
    block = 'clad'
  []  
[]

[AuxVariables]
  [T_wall]
    initial_condition = ${inlet_temperature}
  []
  [heat_source]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = 3e8
    block = 'fuel'
  []
  [q_prime]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = 0.0
    block = 'fuel'
  []
[]

[BCs]
  [cladding_outer_bc]
    type = MatchedValueBC
    variable = T
    v = T_wall
    boundary = '7'
  []
[]

[Executioner]
  type = Transient
[]

[AuxKernels]
  [q_prime]
    type = SpatialUserObjectAux
    variable = q_prime
    user_object = q_prime_uo
    block = 'fuel'
  []
[]

[UserObjects]
  [q_prime_uo]
    type = NearestPointLayeredIntegral
    variable = heat_source
    block = 'fuel'
    direction = z
    points_file = 'pin_centers.txt'
    num_layers = ${n_layers}
    execute_on = 'transfer timestep_end'
  []
[]


[Postprocessors]
  [conduction_power_integral]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    block = 'fuel'
    execute_on = 'transfer'
  []
[]

[Outputs]
  exodus = true
[]`

`!include common.i

[Mesh]
  [file]
    type = FileMeshGenerator
    file = mesh_in.e
  []
[]

[Problem]
  type = OpenMCCellAverageProblem
  power = ${power}
  lowest_cell_level = 2
  scaling = 100
  debug = true
  temperature_blocks = 'helium fuel clad lead'
  density_blocks = 'lead'
  xml_directory = '/home/cardinal/moose_test/subchannel_coupling'

  [Tallies]
    [power]
      type = CellTally
      score = 'kappa_fission'
    []
  []
[]

[ICs]
  [temp]
    type = ConstantIC
    variable = temp
    value = ${inlet_temperature}
  []
  [density]
    type = ConstantIC
    variable = density
    value = ${fparse (11441 - 1.2795 * inlet_temperature)}
  []
[]

[AuxVariables]
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []
  [cell_density]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []
  [cell_density]
    type = CellDensityAux
    variable = cell_density
  []
[]

[Executioner]
  type = Transient
  num_steps = 3
[]

[Outputs]
  exodus = true
[]

[MultiApps]
  [conduction]
    type = TransientMultiApp
    input_files = 'conduction.i'
    sub_cycling = true
    execute_on = timestep_end
  []
  [subchannel]
    type = FullSolveMultiApp
    input_files = 'subchannel.i'
    max_procs_per_app = 1
    execute_on = timestep_begin
  []
[]

[Transfers]
  [power_to_conduction]
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = conduction
    source_variable = kappa_fission
    variable = heat_source
    from_postprocessors_to_be_preserved = openmc_power_integral
    to_postprocessors_to_be_preserved = conduction_power_integral
  []
  [solid_temperature_from_conduction]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = conduction
    source_variable = T
    variable = temp
    to_blocks = 'fuel clad helium'
  []

  [linear_heat_rate_to_subchannel]
    type = MultiAppGeneralFieldNearestLocationTransfer
    variable = q_prime
    source_variable = q_prime
    from_multi_app = conduction
    to_multi_app = subchannel
    greedy_search = true
    use_bounding_boxes = false
    to_blocks = 'fuel_pins'
  []
  [fluid_temperature_from_subchannel]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = T
    variable = temp
    from_multi_app = subchannel
    greedy_search = true
    use_bounding_boxes = false
    to_blocks = 'lead'
    from_blocks = 'subchannel'
  []
  [fluid_density_from_subchannel]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = rho
    variable = density
    from_multi_app = subchannel
    greedy_search = true
    use_bounding_boxes = false
    to_blocks = 'lead'
    from_blocks = 'subchannel'
  []
  [clad_surface_temperature_to_conduction]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = subchannel
    to_multi_app = conduction
    source_variable = Tpin
    variable = T_wall
  []
[]

[Postprocessors]
  [openmc_power_integral]
    type = ElementIntegralVariablePostprocessor
    variable = kappa_fission
    execute_on = 'transfer timestep_end'
  []
[]`

`!include common.i

temperature = ${fparse 0.5 * (outlet_temperature + inlet_temperature)}
Cp = ${fparse 145}
total_mdot = ${fparse power / Cp / (outlet_temperature - inlet_temperature)}

flow_area = ${fparse 3*sqrt(3)/2 * (edge_length*1e-2)^2 - 169 * 3.14 * (clad_outer_diameter*1e-2)^2/4}
mass_flux_in = ${fparse total_mdot / flow_area}

[TriSubChannelMesh]
  [subchannel]
    type = SCMTriSubChannelMeshGenerator
    nrings = 8
    n_cells = 100
    flat_to_flat = ${fparse edge_length * sqrt(3)* 1e-2}
    heated_length = ${fparse height * 1e-2}
    pin_diameter = ${fparse clad_outer_diameter * 1e-2}
    pitch = ${fparse fuel_pin_pitch * 1e-2}
    dwire = ${fparse wire_diameter * 1e-2}
    hwire = ${fparse wire_pitch * 1e-2}
  []
  [fuel_pins]
    type = SCMTriPinMeshGenerator
    input = subchannel
    nrings = 8
    n_cells = 100
    heated_length = ${fparse height * 1e-2}
    pitch = ${fparse fuel_pin_pitch * 1e-2}
  []
[]

[AuxVariables]
  [mdot]
    block = subchannel
  []
  [SumWij]
    block = subchannel
  []
  [P]
    block = subchannel
    initial_condition = ${P_out}
  []
  [DP]
    block = subchannel
  []
  [h]
    block = subchannel
  []
  [T]
    block = subchannel
    initial_condition = 693
  []
  [rho]
    block = subchannel
  []
  [S]
    block = subchannel
  []
  [w_perim]
    block = subchannel
  []
  [mu]
    block = subchannel
  []
  [Tpin]
    block = fuel_pins
    initial_condition = ${inlet_temperature}
  []
  [q_prime]
    block = fuel_pins
    initial_condition = ${fparse power/169/(height * 1e-2)}
  []
  [Dpin]
    block = fuel_pins
    initial_condition = ${fparse 1e-2 * clad_outer_diameter}
  []
  [displacement]
    block = fuel_pins
  []
[]

[ICs]
  [S_IC]
    type = SCMTriFlowAreaIC
    variable = S
  []
  [w_perim_IC]
    type = SCMTriWettedPerimIC
    variable = w_perim
  []
  [Viscosity_ic]
    type = ViscosityIC
    variable = mu
    p = ${P_out}
    T = T
    fp = lead
  []
  [rho_ic]
    type = RhoFromPressureTemperatureIC
    variable = rho
    p = ${P_out}
    T = T
    fp = lead
  []
  [h_ic]
    type = SpecificEnthalpyFromPressureTemperatureIC
    variable = h
    p = ${P_out}
    T = T
    fp = lead
  []
[]


[FluidProperties]
  [lead]
    type = LeadFluidProperties
  []
[]

[Problem]
  type = TriSubChannel1PhaseProblem
  fp = lead
  n_blocks = 1
  P_out = ${P_out}
  CT = 1.0
  compute_density = true
  compute_viscosity = true
  compute_power = true
  implicit = true
  segregated = false
  staggered_pressure = false
  monolithic_thermal = false
  interpolation_scheme = upwind
  P_tol = 1.0e-5
  T_tol = 1.0e-5
[]

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

[Postprocessors]
  [power]
    type = SCMPinPowerPostprocessor
  []
  [inlet_temp]
    type = SCMPlanarMean
    variable = T
    height = 0.0
  []
  [outlet_temp]
    type = SCMPlanarMean
    variable = T
    height = ${fparse height * 1e-2}
  []
  [Total_power_SCMPinPowerPostprocessor]
    type = SCMPinPowerPostprocessor
  []
  [Total_power_SCMTHPowerPostprocessor]
    type = SCMTHPowerPostprocessor
  []
  [total_mdot]
    type = SCMPlanarMean
    variable = mdot
    height = 0.0
  []
  [mean_S]
    type = SCMPlanarMean
    variable = S
    height = 0.0
  []   
[]

[Executioner]
  type = Steady
[]

[Outputs]
  exodus = true
[]`

[GiudGiud]

Hello

did you check visually that the power distribution is well aligned with the bottom and top of the heated region? Using the exodus output

did you account for the difference in length units?

[GiudGiud]

I think you might be transferring q’’’ (volumetric heat rate) into q’ (linear heat rate)?

the centroids to nodes warning will need to be addressed in a second time

[albottar]

Ok, I found a way how to deal with centroids (with ProjectAux). About your guess with q''' - I think you're right. Just checked kappa-fission from openmc and solid heat_source, they're exactly the same.

q_prime for subchannel pins peaks at 2.8e3, and power (partial), as seen before, is 3.4e5. Just need to find, where exactly mistake is.

[albottar]

Yes, I've made a little research - I'm losing about 95% of heating - SCMPimPowerProcessor (total energy transfered to lead) shows about 3.4e5, by simple math its 95% loss.

UPD: Idk what is going on with my Internet...

[UserObjects]
  [q_prime_uo]
    type = NearestPointLayeredIntegral
    variable = heat_source
    block = 'fuel'
    direction = z
    points_file = 'pin_centers.txt'
    num_layers = ${n_layers}
    execute_on = 'transfer timestep_end'
  []
[]

How can I rewrite it to get linear power?

[GiudGiud]

You can use a SpatialUserObjectAux to get this into an auxvariable, then use a ParsedAux to do the normalization operations

[albottar]

[AuxKernels]
  [q_prime_raw_from_uo]
    type = SpatialUserObjectAux
    variable = q_prime_raw
    user_object = q_prime_uo
    block = 'fuel'
    execute_on = 'timestep_end'
  []
  [q_prime_from_raw]
  type = ParsedAux
  variable = q_prime
  expression = 'q_prime_raw / ${fparse height*1e-2/n_layers}'
  coupled_variables = 'q_prime_raw'
  block = 'fuel'
  execute_on = 'timestep_end'
  []
[]

[UserObjects]
  [q_prime_uo]
    type = NearestPointLayeredIntegral
    variable = heat_source
    block = 'fuel'
    direction = z
    points_file = 'pin_centers.txt'
    num_layers = ${n_layers}
    execute_on = 'timestep_end'
  []
[]

Big handshake from the research reactor staff! Thank you, GiudGiud!