terrapower / armi

An open-source nuclear reactor analysis automation framework that helps design teams increase efficiency and quality
https://terrapower.github.io/armi/
Apache License 2.0
231 stars 89 forks source link

component.density and component.materials.density are not the same #1436

Open keckler opened 1 year ago

keckler commented 1 year ago

Long story short, if you build an assembly with UraniumOxide in it and then expand components of that assembly, you will see that the total mass is different depending on how you calculate it. My belief is that this is caused by the UraniumOxide.density() and UraniumOxide.linearExpansion() and/or UraniumOxide.linearExpansionPercent() methods not being completely aligned with each other. Here is an example of the problem.

Let's build a very simple assembly with just one component, made of UraniumOxide.

import copy

import armi; armi.configure()

from armi.reactor.converters.axialExpansionChanger import AxialExpansionChanger
from armi.reactor import grids, components, assemblies, blocks
from armi.reactor.flags import Flags

bFuel = blocks.HexBlock("fuel", height=1)
bDummy = blocks.HexBlock("dummy", height=1)

bFuel.add(components.Circle("fuel", "UraniumOxide", Tinput=25, Thot=25, od=2, id=0, mult=1))
bDummy.add(components.Hexagon("dummy", "Void", Tinput=25, Thot=25, op=20, ip=0, mult=1))
bDummy.p.flags = Flags.DUMMY
bFuel.setAxialExpTargetComp(bFuel[0])
bDummy.setAxialExpTargetComp(bDummy[0])

a = assemblies.HexAssembly("fuel")
a.spatialGrid = grids.axialUnitGrid(2)
a.spatialGrid.armiObject = a

a.add(copy.deepcopy(bFuel))
a.add(copy.deepcopy(bDummy))

a.calculateZCoords()
a.reestablishBlockOrder()

We can now change the temperature of the component and perform full expansion:

axc = AxialExpansionChanger(detailedAxialExpansion=True)
axc.setAssembly(a)
dT = 500
axc.expansionData.updateComponentTemp(a[0][0], a[0][0].temperatureInC + dT)
axc.expansionData.computeThermalExpansionFactors()
axc.axiallyExpandAssembly()

If we now query the mass of that component in two separate ways, we get different answers:

a[0][0].getMass()
34.319448492499255

a[0][0].getVolume() * a[0][0].material.density(Tc=25+dT)
34.476302875694955

That clearly is not correct. This difference is small, but it presents a problem when one is trying to show that mass is conserved between their blueprints and the actual instantiated ARMI model.

If we repeat the process for UZr instead, we get the same component mass when calculated in the two ways shown:

a[0][0].getMass()
50.34706524380068

a[0][0].getVolume() * a[0][0].material.density(Tc=25+dT)
50.34706524380066

The reason that UZr does this correctly, but UraniumOxide does not, is because UZr defines only the linearExpansionPercent() method and lets density() be calculated relative to that. In UraniumOxide we define both explicitly, and thus they are able to be not perfectly aligned with each other.

I do not know if this problem extends to other material classes, but I imagine that it might. We currently have no guardrails to prevent this sort of thing.

keckler commented 1 year ago

This is probably a specific instance of https://github.com/terrapower/armi/issues/865

john-science commented 2 months ago

Okay, the problem here is clearly that this is not how we calculate mass:

a[0][0].getVolume() * a[0][0].material.density(Tc=25+dT)

We calculate mass starting at the number densities of each nuclide:

https://github.com/terrapower/armi/blob/9dc6760de7da3800cbe98e83066cf5f343e81a6a/armi/reactor/components/component.py#L829-L836

The difference is, of course, that you did this:

c.getVolume() * c.material.density(Tc=T)

when you should have done this:

c.getVolume() * c.density()

When I run this second line, I get the exact result. It matches perfectly. Does that make sense?

keckler commented 2 months ago

Yeah, it has only recently sunk in for me that c.material.density() should pretty much never be used. This has always struck me as not obvious to a user, and thus dangerous.

But I think still, for a component that has experienced no burnup and no physics manipulation of any kind outside of thermal expansion (as the component in my example), wouldn't we expect for the material object's density() to be in agreement with the expanded component? I guess my impression was that we only don't expect that after "physics" has been performed upon the model.

If I'm wrong there, can you explain to me why not?

john-science commented 2 months ago

I'm not sure what you would want to do differently.

People DO use (and need) many other things from a component's material

Those are all real examples taken from the wild. And there are many more. The .material of each Component needs to be accessible.

I believe the general idea was supposed to be that users could do this:

c.density()

but this fails:

c.material.density()

And that should stop people from getting confused. BUT, okay, confusion was had. That's no good. I don't immediately see an improvement that is shorter term than just rewrite the materials library, which will hopefully happen within the next year.

keckler commented 2 months ago

but this fails:

c.material.density()

And that should stop people from getting confused.

Can you clarify what you mean by this? As shown in my example in the initial post above, that type of call does not fail.