Negative and very small values of average fission heat, core end temperature, group1 current, etc.

[Mateusz Pater]

Dear Moltres users,

I noticed that from time to time, when I use different sets of group constants, the computations produce negative and very small numbers as mentioned in the subject. Such simulations crash sometimes, but often they don't, although they are not physical. It is especially the case of MSRE 3d model that is available on github...

Lately I created a new 2d mesh for my own model of a reactor, generated group constants in Serpent, ran an eigenvalue simulation and got a desired keff=1.25. Nevertheless, after setting up proper material values and kernels in a Moltres input file, the simulation couldn't converge. I changed the time step to a smaller one, I also tried out altering the executioner type to PFJNK and NEWTON. From many different alternatives, the one with increased tolerances for iterations succeeded. Nevertheless, after a few time steps, the negative numbers appeared.

Could you please give me a hand with this? I believe you are much more experienced with M&S and you are able to suspect what could be changed in the model to make it work.

Best regards,

Matt

The input file attached:

flow_velocity=23.07

nt_scale=1e13

ini_temp=873

diri_temp=873

[GlobalParams]

  num_groups = 2

  num_precursor_groups = 6

  use_exp_form = false

  group_fluxes = 'group1 group2'

  temperature = temp

  sss2_input = true

  pre_concs = 'pre1 pre2 pre3 pre4 pre5 pre6'

  account_delayed = true

  nt_scale = ${nt_scale}

[]

[Mesh]

  file = 'cmsr.msh'

[]

[Problem]

  coord_type = RZ

[]

[Variables]

  [./temp]

    initial_condition = ${ini_temp}

    scaling = 1e-4

  [../]

[]

[AuxVariables]

  [./power_density]

    order = CONSTANT

    family = MONOMIAL

  [../]

[]

[Precursors]

  [./pres]

    var_name_base = pre

    block = 'fuel'

    outlet_boundaries = 'fuel_tops'

    u_def = 0

    v_def = ${flow_velocity}

    w_def = 0

    nt_exp_form = false

    family = MONOMIAL

    order = CONSTANT

    # jac_test = true

  [../]

[]

[Nt]

  var_name_base = group

  vacuum_boundaries = 'fuel_bottoms fuel_tops moder_bottoms moder_tops outer_wall'

  create_temperature_var = false

  scaling = 1e-4

  pre_blocks = 'fuel'

[]

[Kernels]

  # Temperature

  [./temp_time_derivative]

    type = MatINSTemperatureTimeDerivative

    variable = temp

  [../]

  [./temp_source_fuel]

    type = TransientFissionHeatSource

    variable = temp

    block = 'fuel'

  [../]

 # [./temp_source_mod]

 #   type = GammaHeatSource

 #   variable = temp

 #   block = 'moder'

 #   average_fission_heat = 'average_fission_heat'

 #   gamma = gamma_func

 # [../]

  [./temp_diffusion]

    type = MatDiffusion

    D_name = 'k'

    variable = temp

  [../]

  [./temp_advection_fuel]

    type = ConservativeTemperatureAdvection

    velocity = '0 ${flow_velocity} 0'

    variable = temp

    block = 'fuel'

  [../]

[]

[BCs]

  #[./temp_diri_cg]

  #  boundary = 'fuel_bottoms outer_wall'

  #  type = FlexiblePostprocessorDirichletBC

  #  postprocessor = coreEndTemp

  #  offset = -27.8

  #  variable = temp

  #[../]

    [./temp_diri_cg]

     boundary = 'moder_bottoms fuel_bottoms outer_wall'

     type = FunctionDirichletBC

     function = 'temp_bc_func'

     variable = temp

   [../]

  [./temp_advection_outlet]

    boundary = 'fuel_tops'

    type = TemperatureOutflowBC

    variable = temp

    velocity = '0 ${flow_velocity} 0'

  [../]

[]

[AuxKernels]

  [./fuel]

    block = 'fuel'

    type = FissionHeatSourceTransientAux

    variable = power_density

  [../]

 # [./moderator]

 #   block = 'moder'

 #   type = ModeratorHeatSourceTransientAux

 #   average_fission_heat = 'average_fission_heat'

 #   variable = power_density

 #   gamma = gamma_func

 # [../]

[]

[Functions]

  [./temp_bc_func]

    type = ParsedFunction

    value = '${ini_temp} - (${ini_temp} - ${diri_temp}) * tanh(t/1e-2)'

  [../]

 # [./gamma_func]

 #   type = ParsedFunction

 #   value = '${gamma_frac} * pi^2 / 4 * cos(pi * x / (2. * ${R})) * sin(pi * y / ${H})'

 # [../]

[]

[Materials]

  ........................

[]

[Executioner]

  type = Transient

  end_time = 10000

  nl_rel_tol = 1e-2

  nl_abs_tol = 6e-2

  solve_type = 'PJFNK'

  line_search = none

  petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor'

  petsc_options_iname = '-pc_type'

  petsc_options_value = 'lu'

  nl_max_its = 30

  l_max_its = 100

  dtmin = 1e-5

  [./TimeStepper]

    type = PostprocessorDT

    postprocessor = limit_k

    dt = 1e-3

  [../]

[]

[Preconditioning]

  [./SMP]

    type = SMP

    full = true

  [../]

[]

[Postprocessors]

  [./group1_current]

    type = IntegralNewVariablePostprocessor

    variable = group1

    outputs = 'console csv'

  [../]

  [./group1_old]

    type = IntegralOldVariablePostprocessor

    variable = group1

    outputs = 'console csv'

  [../]

  [./multiplication]

    type = DivisionPostprocessor

    value1 = group1_current

    value2 = group1_old

    outputs = 'console csv'

  [../]

  [./temp_fuel]

    type = ElementAverageValue

    variable = temp

    block = 'fuel'

    outputs = 'csv console'

  [../]

  [./temp_moder]

    type = ElementAverageValue

    variable = temp

    block = 'moder'

    outputs = 'csv console'

  [../]

  [./average_fission_heat]

    type = AverageFissionHeat

    execute_on = 'linear nonlinear'

    outputs = 'csv console'

    block = 'fuel'

  [../]

  [./coreEndTemp]

    type = SideAverageValue

    variable = temp

    boundary = 'fuel_tops'

    outputs = 'csv console'

    execute_on = 'linear nonlinear'

  [../]

  [./limit_k]

    type = LimitK

    execute_on = 'timestep_end'

    k_postprocessor = multiplication

    growth_factor = 1.2

    cutback_factor = .4

    k_threshold = 1.5

  [../]

[]

[Outputs]

  perf_graph = true

  print_linear_residuals = true

  csv = true

  exodus = true

[]

[Debug]

  show_var_residual_norms = true

[]

[Park Sun Myung, -]

Hi Matt,

May I ask how small your time steps are? For some of my models, I've had to start from t = 1e-6 depending on the initial conditions.

If the initial temperatures are too different from the steady state temperatures, the initial flux response tends to be very large due to the strong temperature feedback.

Best,

Sun Myung

--
You received this message because you are subscribed to the Google Groups "moltres-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moltres-user...@googlegroups.com.
To post to this group, send email to moltre...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/moltres-users/50146303-aa1f-4b61-b412-a2d81bb0936c%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Mateusz Pater]

Hi, thank you for your reply.

In the above example I started with dt=1e-3 and dtmin=1e-5 and subsequently I am changing the numbers for lower and lower.

Is there a way to find the steady state temperature? The initial conditions I imposed come from the design of a reactor and are, like in MSRE examples, set as the core inlet temperatures.

To unsubscribe from this group and stop receiving emails from it, send an email to moltre...@googlegroups.com.

[Park Sun Myung, -]

Hi Matt,

Hi Matt,

I think that imposing the inlet temperatures throughout the reactor would result in a large excess reactivity. You could try starting at higher temperatures. As for the steady state temperature, it would have to come from other literature. From my own experience, it is generally more difficult to obtain convergence from using the steady state executioner as opposed to the transient executioner.

Sun Myung

---

From: moltre...@googlegroups.com <moltre...@googlegroups.com> on behalf of Mateusz Pater <mateusz...@gmail.com>
Sent: Tuesday, June 18, 2019 4:27 AM
To: moltres-users
Subject: Re: Negative and very small values of average fission heat, core end temperature, group1 current, etc.

 

To unsubscribe from this group and stop receiving emails from it, send an email to moltres-user...@googlegroups.com.

To post to this group, send email to moltre...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/moltres-users/e36f6553-41a2-4168-80a2-4ef329412523%40googlegroups.com.

[Huff, Kathryn D]

I don't see anything obvious, but I think this kind of behavior may be related to the kernels needing some stabilization scheme for your particular regime, or... potentially you could help to reach the solution within the expected range through judicious use of nt_scaling. 

--
You received this message because you are subscribed to the Google Groups "moltres-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moltres-user...@googlegroups.com.
To post to this group, send email to moltre...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/moltres-users/50146303-aa1f-4b61-b412-a2d81bb0936c%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

Kathryn Huff

Assistant Professor

Nuclear, Plasma, and Radiological Engineering

The University of Illinois at Urbana-Champaign

Email: kdh...@illinois.edu 

Preferred Pronoun: she/her

Individual site: http://kdhuff.npre.illinois.edu/

Research group site: http://arfc.npre.illinois.edu/ 

Office: 118 Talbot Laboratory, 104 S. Wright St. 

Book a meeting with me: https://katyhuff.youcanbook.me/

[Gavin Ridley]

Sounds to me like you may need what is every numerical analyst’s best friend: under-relaxation. In MOOSE terms, you may want to try a damper on the nonlinear iterations.

You mentioned that you needed increased tolerances. That definitely sounds like an issue to be resolved through scaling variables.

As for stabilization of the fluxes, if you’re using a diffusion model of neutronics, stabilization may not help. Almost all stabilization schemes add some sort of false diffusion, but this problem is already pure diffusion.

Lastly, is your temperature feedback coefficient negative in both salt and graphite over all temperatures? If not, you’re not guaranteed to have a steady state.

Gavin Ridley

To view this discussion on the web visit https://groups.google.com/d/msgid/moltres-users/CAFoDED3gYYOnt2a6xGOJTEaSOjJqq4LTNMVDZLVPVAVicpRM1Q%40mail.gmail.com.

[Alexander Lindsay]

We've added an automatic scaling option in MOOSE that scales your residuals and Jacobians in such a way that the most important diagonal terms for each variable are unity.

You can turn this option on using `automatic_scaling=true` in your `Executioner` block.

To view this discussion on the web visit https://groups.google.com/d/msgid/moltres-users/DE4E5A28-AAEC-44F4-B863-CE8E84B4D239%40vols.utk.edu.

[Alexander Lindsay]

Also if your physics is going to change a lot over the course of a transient, you can choose to automatically scale at each time step. Note that this won't currently work with TimeIntegrators that use old solution information, but it works with the default time integrator that you get if you don't specify anything in your input file (ImplicitEuler).

You can scale on every time step by setting `compute_scaling_once=false` in your `Executioner` block.

You'll need a new MOOSE in order to get these capabilities.

[Mateusz Pater]

Thank you all, the problem was apparently in the scaling of variables.

The new MOOSE functionality sounds great - it will solve a lot of issues! Good work!

Best regards,

Matt