Burnup calculation of a 1/8 assembly

[Yue JIN]

Hi all,

I ran a depletion calculation of a pwr assembly. And the results are in good agreement with deterministic code.

I wonder if I can divide an assembly by 8 to reduce the standard deviation of tallies. But it seems that standard deviation doesn't have a noticeable change. Why?

This time, I used a triangular prism instead of a box with reflective boundarys to contain the 17*17 lattice (thus some pins are cut). Fuel.volume is divided by eight while Intergrator.powerdensity remains unchanged.

 

The results become wrong as shown in the figure on the left below, keff calculated by OpenMC decreases faster than deterministic code. However, if I divide both Fuel.volume and Intergrator.powerdensity by 8, the results are more reasonable as shown on the right below (although still wrong). How can this happen while total power is only 1/64 of what it was before?  Does the cutting of some fuel pins affect the evolution of nuclide densities within fuel in these pins? 

 

Regards,

Yue

[Andrew Johnson]

Hi Yue,

I'll address the depletion related questions specifically, as I know handling symmetry and depletion can be tricky.

First, it looks like all three cases give slightly different trends for your multiplication factor. It could be an artifact of plotting / rendering in the browser, but your first model hits k < 1.125 at day 250, while the first pass at 1/8th symmetry reaches slightly greater than 1.12 at t=250 days. Then, in what is presented as the "right" figure out to day 25, the final k value for ROBIN looks to be around 1.23, which appears lower than the two previous cases. This leads me to think there could be some issue in modeling these symmetric models.

Nonetheless, on to OpenMC. There are two key candidates for errors in my perspective: material divisions and volumes. Looking at your first model, it appears that the outer ring of fuel pins is segmented into halves, with quarter pins at the corner. With reflective boundary conditions, this is fine from a transport perspective. Similarly, in the 1/8 model, some pins are (understandably and correctly) cut in half. Nothing wrong from a transport perspective but this can be a source of error for depletion. This is dependent on how you represent the materials as well.

If you model each fuel pin as the same material (suggested by your lattice images), and do not apply the "diff_burnable_mats" option for the "openmc.deplete.Operator", then the total volume of your fuel is correct, provided you determined the total volume of all your fuel pins and applied that to the fuel. However, this will cause the entire fuel to be depleted with an average energy spectrum, while each individual pin sees a truly unique spectrum. This affects the one-group reaction rates necessary for depletion, which in turn will change the overall evolution of you materials. If your deterministic code is depleting the entire fuel as a single material as well, then this is not a source of error.

Issues with material volumes can arise when applying symmetry and using this automated material differentiation. If you have N fuel pins each with the same material, and that material has a volume V, N unique materials will be created, each with a volume V/N. If you have a pin that is bisected in half, however, the true volume of this pin should be V/(2N), and thus the total number of atoms in that material will be off by a factor of 2. This could be a driver behind the increasing difference between ROBIN and OpenMC in the second figure depleted out to 250 days.

Finally, I am surprised to hear that dividing both the fuel volume and power density by 8 gave you better agreement with ROBIN. Is this supposed to be the figure depleted out to 25 days? If you have a full assembly generating a power of W, then (if volumes are properly handled), the 1/8 model would generate a power of 1/8 as well. Yet, the power density should remain the same. The OpenMC line in the figure depleted out to 25 days appears to be greatly over depleted.

Could you clarify what figure corresponds to what model and power density?

Regards,

Andrew

[Yue JIN]

Thank you, Andrew. If I understand you correctly, with “diff_burnable_mats” on, I can’t assign a single “total volume” to fuel material in this case. If there are two cells with different volumes in a geometry, even though the material filled in these two cells is the same, I still need to define two materials. For example, material_1.volume = 10 and material_2.volume =40 ?

Sorry, the two bottom pictures have swapped their positions. I posted several times, but the layout was never correct. I think google should add a “preview” feature.

I’m also surprised that dividing both volume and power density gave better agreement. Maybe I should not divide the total volume by 8? Because I built the second model by directly cutting a 17*17 lattice with a triangular prism. I didn’t build this model by defining a smaller lattice. Whether I cut or not, the volume of each pin is always V/(17*17)?

My last question is why dividing a model by 8 doesn’t reduce the standard deviation. I ran simulation of a single state point (thus no depletion) for two models. The number of particles and batches is the same for two models （200*20000）. But there is no obvious difference between the standard deviations of eigenvalue.

I will retry with a manually constructed 1/8 model instead of directly cutting a complete assembly.

 

Regards,

Yue.

[Andrew Johnson]

Yue,

You are correct in regard to how the "diff_burnable_mats" option would affect building you model. What I have done in the past is create specific materials for pins that represent half, full and quarter pins, with identical isotopics but with adjusted volumes. Then, I can assign these three materials to the correct pins so that the volumes are still properly handled.

Ah, you built a full 17x17 (fuel pins and all) and then cut it with a triangle? This explains the power normalization issue. Let's again say the full assembly produces a power of W with power density D. We ideally want our octant yield a power of W/8 with the same power density. When you provide the power density to OpenMC, it computes the power based on all the materials contained in your problem. What I believe happens is the fuel pins that are not included in the simulation of neutrons (e.g. outside your triangle) contribute to the total heavy metal mass of your octant model, giving you the same heavy metal as in your full assembly model. This in turn, gives your assembly a power of W, rather than the W/8 we want. But we still have to account for the remaining 1/8 scaling, since dividing your power density and fuel volume by 8 (for a total reduction 1/8 * 1/8 = 1/64) yielded closer results to ROBIN.

When OpenMC runs the transport simulation, reaction rates in all burnable materials are tallied, which include those outside your triangle. These will contribute zero power (through the fission q mode or energy deposition), meaning the tallies in your triangle are further increased, since the pins in the simulated eighth must account for the power of a full assembly. You might be able to verify this by looking at the number of burnable materials present in the depletion result file. The full assembly should have around 264 materials (17 x 17 lattice with 25 empty cells), while your octant model should have 33 materials, or 1/8th the original number of fuel pins. You can also use "len(Operator.burnable_mats)" rather than digging in to the output files.

You can still make a 17x17 (or even a square quarter) lattice and place it inside a bounding triangle, which might make modeling easier. I would recommend that you fill the positions that do not exist in the final geometry (fuel and non-fuel pins alike) with a dummy universe, like a pure water pin. Since these pins don't really exist in your model, it doesn't matter (from a neutronics perspective) what you put there. The only differences will be the power normalization (as seen here) and the memory footprint of the simulation.

The volume assigned to your fuel pins should align with the true volume (or cross-sectional area for 2d) of your geometry. If you are modeling 1/8 of your problem, then the total volume of fuel pins should add to 1/8 that of your full assembly. Again, we want to have the same power density, and the octant model produces 1/8 the power, and thus should contain 1/8 the volume of fuel. Distributed carefully across full and half pins, as mentioned earlier.

Finally, I am not completely sure why the standard deviation of k doesn't change, but I have a guess. Since k is computed over your entire problem, cutting the geometry and applying symmetry doesn't change the number of times a neutron contributes to the multiplication factor. If you were to produce a more localized tally, like flux in a specific fuel pin, you would likely have more neutrons contributing to that tally in your symmetric model, and thus the uncertainty would likely be decreased.

Regards,

Andrew

[Yue JIN]

Thank you very much, Andrew! Your explanation of  deviation of k sounds reasonable. I think in the symmetric model I don’t need to set as many inactive batches as in a complete assembly.

Regards,

Yue.