diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst
index e7d6d320f1e..3034826bddd 100644
--- a/docs/source/pythonapi/model.rst
+++ b/docs/source/pythonapi/model.rst
@@ -32,6 +32,7 @@ Composite Surfaces
    openmc.model.RectangularParallelepiped
    openmc.model.RectangularPrism
    openmc.model.RightCircularCylinder
+   openmc.model.Vessel
    openmc.model.XConeOneSided
    openmc.model.YConeOneSided
    openmc.model.ZConeOneSided
diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index df290329647..e9458f70658 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -40,9 +40,11 @@ def boundary_type(self):
     def boundary_type(self, boundary_type):
         # Set boundary type on underlying surfaces, but not for ambiguity plane
         # on one-sided cones
+        classes = (XConeOneSided, YConeOneSided, ZConeOneSided, Vessel)
         for name in self._surface_names:
-            if name != 'plane':
-                getattr(self, name).boundary_type = boundary_type
+            if isinstance(self, classes) and name.startswith('plane'):
+                continue
+            getattr(self, name).boundary_type = boundary_type
 
     def __repr__(self):
         return f"<{type(self).__name__} at 0x{id(self):x}>"
@@ -89,8 +91,8 @@ class CylinderSector(CompositeSurface):
         counterclockwise direction with respect to the first basis axis
         (+y, +z, or +x). Must be greater than :attr:`theta1`.
     center : iterable of float
-       Coordinate for central axes of cylinders in the (y, z), (x, z), or (x, y)
-       basis. Defaults to (0,0).
+        Coordinate for central axes of cylinders in the (y, z), (x, z), or (x, y)
+        basis. Defaults to (0,0).
     axis : {'x', 'y', 'z'}
         Central axis of the cylinders defining the inner and outer surfaces of
         the sector. Defaults to 'z'.
@@ -118,7 +120,7 @@ def __init__(self,
                  r2,
                  theta1,
                  theta2,
-                 center=(0.,0.),
+                 center=(0., 0.),
                  axis='z',
                  **kwargs):
 
@@ -1837,3 +1839,96 @@ def __init__(self, center_base: Sequence[float], axis: Sequence[float],
 
     def __neg__(self) -> openmc.Region:
         return +self.plane_bottom & -self.plane_top & -self.cone
+
+
+class Vessel(CompositeSurface):
+    """Vessel composed of cylinder with semi-ellipsoid top and bottom.
+
+    This composite surface is represented by a finite cylinder with ellipsoidal
+    top and bottom surfaces. This surface is equivalent to the 'vesesl' surface
+    in Serpent.
+
+    .. versionadded:: 0.15.1
+
+    Parameters
+    ----------
+    r : float
+        Radius of vessel.
+    p1 : float
+        Minimum coordinate for cylindrical part of vessel.
+    p2 : float
+        Maximum coordinate for cylindrical part of vessel.
+    h1 : float
+        Height of bottom ellipsoidal part of vessel.
+    h2 : float
+        Height of top ellipsoidal part of vessel.
+    center : 2-tuple of float
+        Coordinate for central axis of the cylinder in the (y, z), (x, z), or
+        (x, y) basis. Defaults to (0,0).
+    axis : {'x', 'y', 'z'}
+        Central axis of the cylinder.
+
+    """
+
+    _surface_names = ('cyl', 'plane_bottom', 'plane_top', 'bottom', 'top')
+
+    def __init__(self, r: float, p1: float, p2: float, h1: float, h2: float,
+                 center: Sequence[float] = (0., 0.), axis: str = 'z', **kwargs):
+        if p1 >= p2:
+            raise ValueError('p1 must be less than p2')
+        check_value('axis', axis, {'x', 'y', 'z'})
+
+        c1, c2 = center
+        cyl_class = getattr(openmc, f'{axis.upper()}Cylinder')
+        plane_class = getattr(openmc, f'{axis.upper()}Plane')
+        self.cyl = cyl_class(c1, c2, r, **kwargs)
+        self.plane_bottom = plane_class(p1)
+        self.plane_top = plane_class(p2)
+
+        # General equation for an ellipsoid:
+        #   (x-x₀)²/r² + (y-y₀)²/r² + (z-z₀)²/h² = 1
+        #   (x-x₀)² + (y-y₀)² + (z-z₀)²s² = r²
+        # Let s = r/h:
+        #   (x² - 2x₀x + x₀²) + (y² - 2y₀y + y₀²) + (z² - 2z₀z + z₀²)s² = r²
+        #   x² + y² + s²z² - 2x₀x - 2y₀y - 2s²z₀z + (x₀² + y₀² + z₀²s² - r²) = 0
+
+        sb = (r/h1)
+        st = (r/h2)
+        kwargs['a'] = kwargs['b'] = kwargs['c'] = 1.0
+        kwargs_bottom = kwargs
+        kwargs_top = kwargs.copy()
+
+        sb2 = sb*sb
+        st2 = st*st
+        kwargs_bottom['k'] = c1*c1 + c2*c2 + p1*p1*sb2 - r*r
+        kwargs_top['k'] = c1*c1 + c2*c2 + p2*p2*st2 - r*r
+
+        if axis == 'x':
+            kwargs_bottom['a'] *= sb2
+            kwargs_top['a'] *= st2
+            kwargs_bottom['g'] = -2*p1*sb2
+            kwargs_top['g'] = -2*p2*st2
+            kwargs_top['h'] = kwargs_bottom['h'] = -2*c1
+            kwargs_top['j'] = kwargs_bottom['j'] = -2*c2
+        elif axis == 'y':
+            kwargs_bottom['b'] *= sb2
+            kwargs_top['b'] *= st2
+            kwargs_top['g'] = kwargs_bottom['g'] = -2*c1
+            kwargs_bottom['h'] = -2*p1*sb2
+            kwargs_top['h'] = -2*p2*st2
+            kwargs_top['j'] = kwargs_bottom['j'] = -2*c2
+        elif axis == 'z':
+            kwargs_bottom['c'] *= sb2
+            kwargs_top['c'] *= st2
+            kwargs_top['g'] = kwargs_bottom['g'] = -2*c1
+            kwargs_top['h'] = kwargs_bottom['h'] = -2*c2
+            kwargs_bottom['j'] = -2*p1*sb2
+            kwargs_top['j'] = -2*p2*st2
+
+        self.bottom = openmc.Quadric(**kwargs_bottom)
+        self.top = openmc.Quadric(**kwargs_top)
+
+    def __neg__(self):
+        return ((-self.cyl & +self.plane_bottom & -self.plane_top) |
+                (-self.bottom & -self.plane_bottom) |
+                (-self.top & +self.plane_top))
diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py
index 963bbe00d19..5504b1c52d5 100644
--- a/tests/unit_tests/test_surface_composite.py
+++ b/tests/unit_tests/test_surface_composite.py
@@ -596,3 +596,52 @@ def test_conical_frustum():
     # Denegenerate case with r1 = r2
     s = openmc.model.ConicalFrustum(center_base, axis, r1, r1)
     assert (1., 1., -0.01) in -s
+
+
+def test_vessel():
+    center = (3.0, 2.0)
+    r = 1.0
+    p1, p2 = -5.0, 5.0
+    h1 = h2 = 1.0
+    s = openmc.model.Vessel(r, p1, p2, h1, h2, center)
+    assert isinstance(s.cyl, openmc.Cylinder)
+    assert isinstance(s.plane_bottom, openmc.Plane)
+    assert isinstance(s.plane_top, openmc.Plane)
+    assert isinstance(s.bottom, openmc.Quadric)
+    assert isinstance(s.top, openmc.Quadric)
+
+    # Make sure boundary condition propagates (but not for planes)
+    s.boundary_type = 'reflective'
+    assert s.boundary_type == 'reflective'
+    assert s.cyl.boundary_type == 'reflective'
+    assert s.bottom.boundary_type == 'reflective'
+    assert s.top.boundary_type == 'reflective'
+    assert s.plane_bottom.boundary_type == 'transmission'
+    assert s.plane_top.boundary_type == 'transmission'
+
+    # Check bounding box
+    ll, ur = (+s).bounding_box
+    assert np.all(np.isinf(ll))
+    assert np.all(np.isinf(ur))
+    ll, ur = (-s).bounding_box
+    assert np.all(np.isinf(ll))
+    assert np.all(np.isinf(ur))
+
+    # __contains__ on associated half-spaces
+    assert (3., 2., 0.) in -s
+    assert (3., 2., -5.0) in -s
+    assert (3., 2., 5.0) in -s
+    assert (3., 2., -5.9) in -s
+    assert (3., 2., 5.9) in -s
+    assert (3., 2., -6.1) not in -s
+    assert (3., 2., 6.1) not in -s
+    assert (4.5, 2., 0.) in +s
+    assert (3., 3.2, 0.) in +s
+    assert (3., 2., 7.) in +s
+
+    # translate method
+    s_t = s.translate((0., 0., 1.))
+    assert (3., 2., 6.1) in -s_t
+
+    # Make sure repr works
+    repr(s)