Hi everyone! I would like to look at making a pebble bed in an annular shape. I’d like to use openmc.model.pack_spheres to place the pebbles in a region bounded by two cylinders. I get the error:
“supported container shapes are rectangular prism, cylinder, sphere, and spherical shell”. I was wondering how difficult it would be to add support for cylindrical shells?
Thanks
I managed to write a cylindrical shell function in triso.py that worked! Havn’t fully error checked it for pushing to main, but here it is for anyone interested:
class _CylindricalShell(_Container):
“”"Cylindrical shell container in which to pack spheres.
Parameters
----------
length : float
Length of the cylindrical container.
radius : float
Outer radius of the cylindrical shell container.
inner_radius : float
Inner radius of the cylindrical shell container.
axis : string
Axis along which the length of the cylinder is aligned.
sphere_radius : float
Radius of spheres to be packed in container.
center : Iterable of float
Cartesian coordinates of the center of the container. Default is
(0., 0., 0.)
Attributes
----------
length : float
Length of the cylindrical container.
radius : float
Outer radius of the cylindrical shell container.
inner_radius : float
Inner radius of the cylindrical shell container.
axis : string
Axis along which the length of the cylinder is aligned.
sphere_radius : float
Radius of spheres to be packed in container.
center : list of float
Cartesian coordinates of the center of the container. Default is
(0., 0., 0.)
shift : list of int
Rolled indices of the x-, y-, and z- coordinates of a sphere so the
configuration is aligned with the correct axis. No shift corresponds to
a cylinder along the z-axis.
cell_length : list of float
Length in x-, y-, and z- directions of each cell in mesh overlaid on
domain.
limits : list of float
Maximum radial distance and minimum and maximum distance in the
direction parallel to the axis where sphere center can be placed.
volume : float
Volume of the container.
"""
def __init__(self, length, radius, inner_radius, axis, sphere_radius, center=(0., 0., 0.)):
super().__init__(sphere_radius, center)
self._shift = None
self.length = length
self.radius = radius
self.inner_radius = inner_radius
self.axis = axis
@property
def length(self):
return self._length
@length.setter
def length(self, length):
self._length = float(length)
self._limits = None
self._cell_length = None
@property
def radius(self):
return self._radius
@radius.setter
def radius(self, radius):
self._radius = float(radius)
self._limits = None
self._cell_length = None
@property
def inner_radius(self):
return self._inner_radius
@inner_radius.setter
def inner_radius(self, inner_radius):
self._inner_radius = float(inner_radius)
self._limits = None
self._cell_length = None
@property
def axis(self):
return self._axis
@axis.setter
def axis(self, axis):
self._axis = axis
self._shift = None
@property
def shift(self):
if self._shift is None:
if self.axis == 'x':
self._shift = [1, 2, 0]
elif self.axis == 'y':
self._shift = [2, 0, 1]
else:
self._shift = [0, 1, 2]
return self._shift
@property
def limits(self):
if self._limits is None:
z0 = self.center[self.shift[2]]
z = self.length/2
sphere_r = self.sphere_radius
r_max = self.radius - sphere_r
if self.inner_radius == 0:
r_min = 0
else:
r_min = self.inner_radius + sphere_r
self._limits = [[z0 - z + sphere_r, r_min], [z0 + z - sphere_r, r_max]] # min and max [z, r] pairs
return self._limits
@limits.setter
def limits(self, limits):
self._limits = limits
@property # leaving this the same as for the cylinder.
def cell_length(self):
if self._cell_length is None:
h = 4*self.sphere_radius
i, j, k = self.shift
self._cell_length = [None]*3
self._cell_length[i] = 2*self.radius/int(2*self.radius/h)
self._cell_length[j] = 2*self.radius/int(2*self.radius/h)
self._cell_length[k] = self.length/int(self.length/h)
return self._cell_length
@property
def volume(self):
return self.length*pi*(self.radius**2 - self.inner_radius**2)
@classmethod
def from_region(self, region, sphere_radius):
check_type('region', region, openmc.Region)
# Assume the simplest case where the cylinder shell volume is the
# intersection of the half-spaces of two cylinders and two planes
if not isinstance(region, openmc.Intersection):
raise ValueError
if any(not isinstance(node, openmc.Halfspace) for node in region):
raise ValueError
if len(region) != 4: # should be two cylinders and two planes
raise ValueError
# Identify the axis that the cylinder lies along
axis = region[0].surface.type[0]
# Make sure the region is composed of a cylinder and two planes on the
# same axis
count = Counter(node.surface.type for node in region)
if count[axis + '-cylinder'] != 2 or count[axis + '-plane'] != 2:
raise ValueError
# Sort the half-spaces by surface type
cyl1, cyl2, p1, p2 = sorted(region, key=lambda x: x.surface.type)
# order the cylinders correctly (cyl1 should be smaller than cyl2)
if cyl1.surface.r > cyl2.surface.r:
cyl1, cyl2 = cyl2, cyl1
# check to make sure cylinders are centered correctly
radius = cyl2.surface.r
inner_radius = cyl1.surface.r
center_check = (cyl1.surface.x0, cyl1.surface.y0, cyl1.surface.z0)
if center_check != (cyl2.surface.x0, cyl2.surface.y0, cyl2.surface.z0):
raise ValueError
# Calculate the parameters for a cylinder along the x-axis
if axis == 'x':
if p1.surface.x0 > p2.surface.x0:
p1, p2 = p2, p1
length = p2.surface.x0 - p1.surface.x0
center = (p1.surface.x0 + length/2, cyl1.surface.y0, cyl1.surface.z0) # TODO cyl1 over cyl2 chosen for no reason
# Calculate the parameters for a cylinder along the y-axis
elif axis == 'y':
if p1.surface.y0 > p2.surface.y0:
p1, p2 = p2, p1
length = p2.surface.y0 - p1.surface.y0
center = (cyl1.surface.x0, p1.surface.y0 + length/2, cyl1.surface.z0) # TODO cyl1 over cyl2 chosen for no reason
# Calculate the parameters for a cylinder along the z-axis
else:
if p1.surface.z0 > p2.surface.z0:
p1, p2 = p2, p1
length = p2.surface.z0 - p1.surface.z0
center = (cyl1.surface.x0, cyl1.surface.y0, p1.surface.z0 + length/2)
# Make sure the half-spaces are on the correct side of the surfaces
if cyl1.side != '+' or cyl2.side != '-' or p1.side != '+' or p2.side != '-':
raise ValueError
# the region is a cylindrical shell, so create container
print('Wow you just made a cylinderical shell!')
return _CylindricalShell(length, radius, inner_radius, axis, sphere_radius, center)
def random_point(self):
ll, ul = self.limits
r = sqrt(uniform(ll[1]**2, ul[1]**2)) # modified to start at the lower radius
t = uniform(0, 2*pi)
i, j, k = self.shift
p = [None]*3
p[i] = r*cos(t) + self.center[i]
p[j] = r*sin(t) + self.center[j]
p[k] = uniform(ll[0], ul[0])
return p
def repel_spheres(self, p, q, d, d_new):
# Moving each sphere distance 's' away from the other along the line
# joining the sphere centers will ensure their final distance is
# equal to the outer diameter
s = (d_new - d)/2
v = (p - q)/d
p += s*v
q -= s*v
# Enforce the rigid boundary by moving each sphere back along the
# surface normal until it is completely within the container if it
# overlaps the surface
ll, ul = self.limits
c = self.center
i, j, k = self.shift
r = sqrt((p[i] - c[i])**2 + (p[j] - c[j])**2)
if r > ul[1]:
p[i] = (p[i] - c[i])*ul[1]/r + c[i] # cartesian "x and y" (after shift)
p[j] = (p[j] - c[j])*ul[1]/r + c[j]
elif r < ll[1]:
p[i] = (p[i] - c[i])*ll[1]/r + c[i] # new section for inner cylinder cartesian "x and y" (after shift)
p[j] = (p[j] - c[j])*ll[1]/r + c[j]
p[k] = np.clip(p[k], ll[0], ul[0])
r = sqrt((q[i] - c[i])**2 + (q[j] - c[j])**2)
if r > ul[1]:
q[i] = (q[i] - c[i])*ul[1]/r + c[i]
q[j] = (q[j] - c[j])*ul[1]/r + c[j]
elif r < ll[1]:
q[i] = (q[i] - c[i])*ll[1]/r + c[i]
q[j] = (q[j] - c[j])*ll[1]/r + c[j]
q[k] = np.clip(q[k], ll[0], ul[0])