Transient simulation

[js-jixu]

Hi, experts.

I want to do a transient nuclear-thermal simulation of a system (single rod for simplicity). I know that we need to use steady-state nuclear thermal simulation to get a solution, which includes neutron flux, precursors, eigenvalue, etc. This steady-state solution is then given to the transient simulation as an initial condition.

When I finished the above work and performed the transient simulation, I found that the neutron flux and powernorm have been decreasing. The specific values can be seen in transient-th-nts.csv, and the initial value (steady-state solution) can be seen in steady-th- nts.csv.

I want to ask is this situation correct? The initial value of powernorm I got in the steady state simulation is 21230, why does powernorm drop so much?

steady2transient.zip

[js-jixu]

I divide nuclear-thermal coupling into nuclear and thermal. The global temperature of the neutron physics field is 600K. First calculate keff for the steady state neutron field, then divide the NSF in xs by keff. This gives us a set of cross-sections that allow the system to reach criticality. I did steady state calculation with modified xs and got keff equal to 1.00000.

Then I used the neutron flux distribution and pres obtained from the steady-state calculation as the initial conditions for the transient calculation, and the modified xs was used for the transient calculation. I set eigenvalue_scaling = 1. It stands to reason that the neutron flux distribution should not change in the transient calculation, because keff is exactly equal to 1, and the whole system is in a critical state. But the results show that the neutron flux and total fission number have been declining. Why is this so?

The attached files have transient and steady state input and output files, modified xs.

modify_xs.zip

[smpark7]

The eigenvalue_scaling parameter scales the neutron source term (fission + delayed) in the neutron diffusion equation. When Moltres runs an eigenvalue calculation with power normalization (e.g. your steady-nts.i file), it normalizes the neutron flux at the end of the simulation to match your desired power output. However, because of how Moltres handles eigenvalue calculations, Moltres does not normalize the delayed neutron precursor distribution. This is why I have the scale_factor = 7.61666e+17 parameter in my cnrs-benchmark phase-2 input files for the precursor variables. Moltres will report a scaling factor at the end of your steady-nts.i simulation. Add this to your pre SolutionFunction in your transient-nts.i file. Let me know if this resolves your issue!

Also, you don't need to import the group fluxes in SolutionUserObject because [Nt] already imports it with init_nts_from_file = true. We use SolutionUserObject instead of init_from_file for [Precursors] because of the normalization issue described above.

[js-jixu]

Thank you, Sun Myung, your suggestion works fine. The transient simulation now works as expected.

Is there any way to increase the precision of scale_factor or output this parameter. This parameter seems to be very important for simulation. I'm concerned that six significant figures may not be precise enough. I have a postprecessor named power_sum to integrate the auxvariable power:

[AuxVariables]
  [power]
    family = MONOMIAL
    order = FIRST
    block = 'fuel'
  []
[]

[AuxKernels]
  [power_auxk]
    type = FissionHeatSourceAux
    variable = power
    tot_fission_heat = powernorm
    power = 21230
  []  
[]

[Postprocessors]
  [power_sum]
    type = ElementIntegralVariablePostprocessor
    variable = power
    block = 'fuel'
    execute_on = linear
  []
[]

But there is a slight difference between the values of power_sum and powernorm in the transient simulation. So I think it may be because the input scale_factor is not accurate enough:

Postprocessor Values:
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+
| time           | bnorm          | group2diff     | group2max      | group2norm     | power_sum      | powernorm      | tot_fissions   |
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+
|   0.000000e+00 |   0.000000e+00 |   0.000000e+00 |   0.000000e+00 |   0.000000e+00 |   0.000000e+00 |   0.000000e+00 |   0.000000e+00 |
|   1.000000e+00 |   1.956703e+15 |   6.925954e+10 |   1.458850e+15 |   1.147499e+17 |   2.123000e+04 |   2.123013e+04 |   6.671718e+14 |
|   2.000000e+00 |   1.956705e+15 |   1.351257e+10 |   1.458851e+15 |   1.147501e+17 |   2.123000e+04 |   2.123015e+04 |   6.671726e+14 |
|   3.000000e+00 |   1.956707e+15 |   1.136042e+10 |   1.458853e+15 |   1.147502e+17 |   2.123000e+04 |   2.123017e+04 |   6.671733e+14 |
|   4.000000e+00 |   1.956709e+15 |   1.001662e+10 |   1.458854e+15 |   1.147503e+17 |   2.123000e+04 |   2.123019e+04 |   6.671738e+14 |
|   5.000000e+00 |   1.956710e+15 |   9.152737e+09 |   1.458855e+15 |   1.147504e+17 |   2.123000e+04 |   2.123021e+04 |   6.671744e+14 |
|   6.000000e+00 |   1.956712e+15 |   8.537862e+09 |   1.458856e+15 |   1.147505e+17 |   2.123000e+04 |   2.123022e+04 |   6.671749e+14 |
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+

[smpark7]

I think the difference you see here is due to your power variable being MONOMIAL. There may be some discretization error from converting power derived from your lagrange-based group fluxes to the discontinuous monomial-based power. Can you try with LAGRANGE?

As for the scale_factor, I think you can calculate it by taking powernorm/bnorm from your steady-nts.i simulation in the final iteration just before normalization.

[js-jixu]

I set power as LAGRANGE variable, but it seems that FissionHeatSourceAux does not support nodal auxiliary variable. The terminal reports that:

*** ERROR ***
The following error occurred in the object "power_auxk", of type "FissionHeatSourceAux".

Nodal AuxKernel 'power_auxk' attempted to reference material property 'fissxs'
Consider using an elemental auxiliary variable for 'power'.

The possible reason is that FissionHeatSourceAux uses computeValue().

For the scale_factor, I used a more precise number. But as you say, it's not the point. powernorm is still a little bigger than the preset power (21230).

[smpark7]

Ah yes I forgot FissionHeatSourceAux only supports elemental variables.

Then, I think power_sum is exactly 21230 because FissionHeatSourceAux calculates the power and normalizes it by powernorm. Therefore, the differences are canceled out.

Can you run the simulation longer to see whether powernorm continues to increase? I think it will increase indefinitely because you set the temperature to be fixed at 600K. We can't get a perfectly stable system because we do not have infinitely precise eigenvalue and scale factors. If you set the temperature at the boundary to 600K and let the body temperature freely vary, I think you can get a stable system from the slight temperature increase => negative feedback. You'll need to define temperature time derivative and diffusion kernels for this approach.

[js-jixu]

Then, I think power_sum is exactly 21230 because FissionHeatSourceAux calculates the power and normalizes it by powernorm. Therefore, the differences are canceled out.

This is right! I think the problem may not be in the error caused by the conversion between the LAGRANGE variable and the MONOMIAL variable. If the problem lies in the error caused by the conversion between the two variables, then the error should be reduced by refining the mesh. I used two meshes, one with three times as many meshes as the other. However, the error didn't decrease.

Then after your inspiration, I think it may be because the eig_check_tol I took in the steady state calculation is not small enough and the scale_factor is not accurate enough. So I compared the results of eig_check_tol=1e-7 and eig_check_tol=1e-10, it shows that the error of eig_check_tol=1e-10 is much smaller than that of eig_check_tol=1e-7.

I attached the data files of eig_check_tol=1e-10 and eig_check_tol=1e-7.

data.zip