diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 1e4b4aa9b..50ad70680 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -387,7 +387,7 @@ class _CylindricalDomain(_Domain): @property def volume(self): - return self.length * pi * self.radius**2 + return self.length*pi*self.radius**2 @length.setter def length(self, length): @@ -493,7 +493,7 @@ class _SphericalDomain(_Domain): @property def volume(self): - return 4/3 * pi * self.radius**3 + return 4/3*pi*self.radius**3 @radius.setter def radius(self, radius): @@ -535,6 +535,125 @@ class _SphericalDomain(_Domain): q *= r_max/r +class _SphericalShellDomain(_Domain): + """Spherical shell container in which to pack particles. + + Parameters + ---------- + radius : float + Outer radius of the spherical shell container. + inner_radius : float + Inner radius of the spherical shell container. + center : Iterable of float + Cartesian coordinates of the center of the container. Default is + [0., 0., 0.] + + Attributes + ---------- + radius : float + Outer radius of the spherical shell container. + inner_radius : float + Inner radius of the spherical shell container. + particle_radius : float + Radius of particles to be packed in container. + center : list of float + Cartesian coordinates of the center of the container. Default is + [0., 0., 0.] + 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 radial distance where particle + center can be placed. + volume : float + Volume of the container. + + """ + + def __init__(self, radius, inner_radius, particle_radius, + center=[0., 0., 0.]): + super().__init__(particle_radius, center) + self.radius = radius + self.inner_radius = inner_radius + + @property + def radius(self): + return self._radius + + @property + def inner_radius(self): + return self._inner_radius + + @property + def limits(self): + if self._limits is None: + self._limits = [self.radius - self.particle_radius, + self.inner_radius + self.particle_radius] + return self._limits + + @property + def cell_length(self): + if self._cell_length is None: + mesh_length = 3*[2*self.radius] + self._cell_length = [x/int(x/(4*self.particle_radius)) + for x in mesh_length] + return self._cell_length + + @property + def volume(self): + return 4/3*pi*(self.radius**3 - self.inner_radius**3) + + @radius.setter + def radius(self, radius): + self._radius = float(radius) + self._limits = None + self._cell_length = None + + @inner_radius.setter + def inner_radius(self, inner_radius): + self._inner_radius = float(inner_radius) + self._limits = None + + @limits.setter + def limits(self, limits): + self._limits = limits + + def random_point(self): + r_max = self.limits[0] + r_min = self.limits[1] + x = (gauss(0, 1), gauss(0, 1), gauss(0, 1)) + r = (uniform(r_min**3, r_max**3)**(1/3)/sqrt(x[0]**2 + x[1]**2 + x[2]**2)) + return [r*s for s in x] + + def repel_particles(self, p, q, d, d_new): + # Moving each particle distance 's' away from the other along the line + # joining the particle 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 particle back along the + # surface normal until it is completely within the container if it + # overlaps the surface + r_max = self.limits[0] + r_min = self.limits[1] + + r = sqrt(p[0]**2 + p[1]**2 + p[2]**2) + if r > r_max: + p *= r_max/r + elif r < r_min: + p *= r_min/r + + r = sqrt(q[0]**2 + q[1]**2 + q[2]**2) + if r > r_max: + q *= r_max/r + elif r < r_min: + q *= r_min/r + + def create_triso_lattice(trisos, lower_left, pitch, shape, background): """Create a lattice containing TRISO particles for optimized tracking. @@ -937,9 +1056,10 @@ def _close_random_pack(domain, particles, contraction_rate): def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, - domain_radius=None, domain_center=[0., 0., 0.], - n_particles=None, packing_fraction=None, - initial_packing_fraction=0.3, contraction_rate=1.e-3, seed=1): + domain_radius=None, domain_inner_radius=None, + domain_center=[0., 0., 0.], n_particles=None, + packing_fraction=None, initial_packing_fraction=0.3, + contraction_rate=1.e-3, seed=1): """Generate a random, non-overlapping configuration of TRISO particles within a container. @@ -955,6 +1075,8 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, Length of the container (if cube or cylinder). domain_radius : float Radius of the container (if cylinder or sphere). + domain_inner_radius : float + Inner radius of the container (if spherical shell). domain_center : Iterable of float Cartesian coordinates of the center of the container. n_particles : int @@ -1018,16 +1140,20 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, """ # Check for valid container geometry and dimensions - if domain_shape not in ['cube', 'cylinder', 'sphere']: + if domain_shape not in ['cube', 'cylinder', 'sphere', 'spherical shell']: raise ValueError('Unable to set domain_shape to "{}". Only "cube", ' - '"cylinder", and "sphere" are ' + '"cylinder", "sphere", and "spherical shell" are ' 'supported."'.format(domain_shape)) if not domain_length and domain_shape in ['cube', 'cylinder']: raise ValueError('"domain_length" must be specified for {} domain ' 'geometry '.format(domain_shape)) - if not domain_radius and domain_shape in ['cylinder', 'sphere']: + if not domain_radius and domain_shape in ['cylinder', 'sphere', + 'spherical shell']: raise ValueError('"domain_radius" must be specified for {} domain ' 'geometry '.format(domain_shape)) + if not domain_inner_radius and domain_shape in ['spherical shell']: + raise ValueError('"domain_inner_radius" must be specified for {} domain ' + 'geometry '.format(domain_shape)) if domain_shape == 'cube': domain = _CubicDomain(length=domain_length, particle_radius=radius, @@ -1038,6 +1164,11 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, elif domain_shape == 'sphere': domain = _SphericalDomain(radius=domain_radius, particle_radius=radius, center=domain_center) + elif domain_shape == 'spherical shell': + domain = _SphericalShellDomain(radius=domain_radius, + inner_radius=domain_inner_radius, + particle_radius=radius, + center=domain_center) # Calculate the packing fraction if the number of particles is specified; # otherwise, calculate the number of particles from the packing fraction.