{'fraction_non_alphanumeric': 0.03447292363857621, 'fraction_numerical': 0.014298752662001825, 'mean_word_length': 4.313497676298242, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mathematics of Floating 3D Printed Objects

We explore the stability of floating objects through mathematical modeling and experimentation. Our models are based on standard ideas of center of gravity, center of buoyancy, and Archimedes' Principle. We investigate a variety of floating shapes with two-dimensional cross sections and identify analytically and/or computationally a potential energy landscape that helps identify stable and unstable floating orientations. We compare our analyses and computations to experiments on floating objects designed and created through 3D printing. In addition to our results, we provide code for testing the floating configurations for new shapes, as well as giving details of the methods for 3D printing the objects. The paper includes conjectures and open problems for further study.


Introduction
Interest in the dynamics of icebergs has been driven by the desire to describe various natural phenomena (including their rolling) as well as practical considerations associated with shipping and protection of offshore structures in arctic environments. Allaire's [Al] study of iceberg stability was motivated by the need to assess mitigation strategies, such as towing of icebergs, to reduce threats posed by icebergs to offshore structures. Allaire identified readily-identifiable above-water characteristics, such as the ratio of waterline width to above-water height, to estimate stability for a menagerie of iceberg shapes -Blocky, Drydock, Dome, Pinnacled, Tabular, Growler. Bailey [B], motivated by similar concerns, examined stability of icebergs in terms of rolling frequency. The potential for iceberg stability considerations to be dynamic even in calm water due to underwater melting/dissolution was considered by Deriabyn & Hjorth [DH]. They were also interested to identify what practical, above-water, observations could be made to predict stability changes driven by underwater changes in iceberg morphology.
Ship design and the design of other man-made floating objects has no doubt driven much scientific and technical work in this field. We make no attempt to review this literature but interested readers may find resources in the work of Mégel & Kliava [MK] and Wilczynski & Diehl [WD]. Historically speaking, scientific thought on this dates back to . See for example Rorres [Ro], a fascinating article on the original work of Archimedes and extensions thereof.
Whereas icebergs and ships are complex three dimensional floating objects, the floating objects that are the focus of the present work are those whose configurations can effectively be characterized in two-dimensions. Specifically, we shall consider 'long' objects whose cross sections are constant for the full length of the object. Such shapes have been the focus of popular online Apps, such as Iceberger [Ice1] and the remixed version [Ice2]. These effective two-dimensional floating objects have proven to be mathematically tractable yet rich in observable phenomena. Prediction of the stable orientations for a long beam with square-cross section and uniform density, for example, was considered by Reid [Re]. Depending on the ratio of the density of the object to the density of the fluid any rotation of the square can be a stable floating orientation (i.e. any orientation from flat side up to corner up). This configuration has been recently revisited both experimentally and analytically by Feigel and Fuzailov [FF]. These authors validated experimentally that for a small range of density ratios near 0.25 (and also similarly near 0.75) the full range of stable orientations can be realized. Another view of this which we discuss in more detail in our work is that owing to the four-fold symmetry of the square, there can be either four or eight stable orientations of the floating square depending on the density ratio. We explore how these results change for what we denote as 'off-center' squares.
The case of a long 'floating plank' of rectangular cross section was the focus of work by Delbourgo [D]. Here, in addition to the density ratio, another parameter -the aspect ratio of the rectangle -appears. Delbourgo identified in this density ratio vs. aspect ratio space the existence of six characteristic floating configurations in terms of (1) long side up or short side up, (2) top side parallel or not parallel to the waterline, and (2) the number of submerged vertices. Further studies have explored other cross sectional shapes and investigated details of the breaking of left-right symmetries of the floating shapes as the density ratio is varied (e.g. see Erdös,Schibler,& Herndon [E1] for square and equilateral triangle cross sections and Erdös, Schibler, & Herndon [E2] for the three-dimensional shapes of the cube, octahedron, and decahedron). Work in this area has focused on homogeneous objects with uniform density. We shall relax this assumption in our present work to include some special classes of non-uniform density in which the center of gravity of the square no longer resides at its centroid (cf. Definition 2.1).
An excellent review of many important mathematical ideas -Archimedes' Principle, Center of Gravity, Center of Buoyancy, and the notion of Metacenter -related to floating objects is the work of Gilbert [G]. One particularly useful concept is that of the potential energy of floating objects. Gilbert shows that if one can compute the potential energy landscape as a function of all possible orientations of the object one can identify stable floating configurations by the locations of local minima of the potential energy function. This is one of the main objectives of our computations as we then use this to identify stable orientations.
The present work shares some of the spirit of the paper by Feigel and Fuzailov [FF] to revisit these questions both theoretically and experimentally. Our experiments, however, are conducted with 3D printed shapes. While some of our theoretical effort has been on objects with square cross section, our methodology is motivated by the recognition that 3D printing offers the opportunity to float objects whose cross sections are effectively limited only by ones own creativity in defining new shapes. An example shape we analyze is the 'Mason M' shown in Figure 1.
The new contributions we make here are (1) predicting floating configurations of an object with square cross section when the center of gravity is not at the object's centroid, (2) predicting floating configurations for objects of general twodimensional polygonal cross section, and (3) verifying the predictions of (1) and (2) as well as classical ones for objects of square cross section with experiments conducted with 3D printed objects. This paper is organized as follows. In Section 2, we introduce definitions and terminology, outline prior results for long floating objects of uniform density with square cross section, and state our analytical results for the case of long floating objects with square cross section with non-uniform density. In Section 3 we give a detailed description of how to 3D print floating objects. This section contains sufficient details and accompanying code so that interested readers would be able to print their own objects and perform their own experiments. Note that we have additionally provided codes and files in a GitHub repository [GIT]. In Section 4 we describe how we obtain experimental results on our 3D printed floating objects. In Section 5 we describe our code for computing stable floating configurations for objects with general polygonal cross sections. In Section 6 we describe the results of our floating experiments and relate them to our theoretical results for several cases involving square cross sections and a selected case of a nontrivial cross section. In Section 7 we give our conclusions and discuss some open problems.

Mathematical Models for Floating Objects
2.1. Definitions. In this section, we give some basic definitions and concepts that will be needed for the remainder of our analysis.
Definition 2.1 (Center of Gravity). If the mass distribution is given by a continuous density function of ρ(x, y, z) within a domain Ω, then the center of gravity can be obtained by (G x where M obj is the object's mass. In the case of uniform density, i.e. ρ is a constant independent of x, y and z, the center of gravity is called the centroid. For a long object of length L with uniform cross section, with uniform density in the long direction, i.e. ρ(x, y, z) = ρ(x, y) is independent of z, the center of gravity is given by with G z = L/2 (relative to one end of the object).
Lemma 2.2. For an object of length L with a polygonal cross-section with uniform constant density ρ, we can compute the area and center of gravity as sums involving only the vertices of the polygon. In particular, let be the vertices of the cross section polygon, oriented counterclockwise. Then the mass of the object is where the area A of the polygon of cross section is given by a result also known as the shoelace formula, and the center of gravity G = (G x , G y ) is given by Proof. According to Green's theorem We can parametrize the line segment between (x k , y k ) and (x k+1 , y k+1 ) by (x, y) = (1 − t)(x k , y k ) + t(x k+1 , y k+1 ) where t ∈ (0, 1). For this line segment, we get dx = (x k+1 − x k )dt and dy = (y k+1 − y k )dt. Now we evaluate the integral from Green's theorem: Now we use this to calculate the full area: Now we turn our attention to the calculation of G x . We use Green's theorem but this time with (f, g) = (0, x 2 /2) in order to satisfy dg dx − df dy = x. Again for the line segment from (x k , y k ) to (x k+1 , y k+1 ) we get Since ρ is constant, we get that G x is equal to ρL/M obj times this integral, but M obj = ρLA, and this gives the factor of 1/A written in the formula above. In a similar way, using Green's theorem with (f, g) = (−y 2 /2, 0), we get This gives the equation for G y stated above.
The proof above relies on Green's theorem, but there are other ways to derive this formula, such as dividing the object into triangles or trapezoids, and combining the corresponding triangle or trapezoid area and center of gravity formulas. Definition 2.3 (Buoyancy). Buoyancy is a force exerted on an object that is wholly or partially submerged in a fluid. The magnitude of this force is equal to the weight of the displaced fluid. Buoyancy relates to the density of the fluid, the volume of the displaced fluid, and the gravitational field; it is independent of the mass and density of the immersed object. The buoyancy force acts vertically upward at the centroid of the displaced volume. The center of buoyancy is given by where V sub is the submerged volume of the object, and Ω sub is the submerged domain. Like in the case of center of gravity, in the uniform cross section case B z = L/2. We use the notation B = (B x , B y ). Note that B is the centroid of the submerged domain.
We now state one of the major results necessary for understanding floating objects.
Theorem 2.4 (Archimedes' Principle). The upward buoyant force exerted on an object wholly or partially submerged is equal to the weight of the displaced fluid. In the absence of other forces, such as surface tension, this can be expressed in the force balance as where g is acceleration due to gravity, ρ f is the density of the fluid, and V sub is the submerged volume of the object.
Observe that (2.4) represents a balanced (net zero) equation of competing forces with the left terms representing gravitational force and right term representing the opposing force of buoyancy. Note that if the object has uniform density, then the mass of the object can be written as M obj = ρ obj V obj . In this case, it follows that For our purposes, Archimedes' Principle determines the appropriate waterline intersections defining a submerged volume whose value relative to the total volume matches the appropriate density ratio. However, it is important to note that satisfying Archimedes' Principle is not a sufficient condition for determining a stable equilibrium. An equilibrium orientation of a floating body occurs when the center of gravity and the center of buoyancy are vertically aligned. If G lies directly below B the equilibrium is stable, whereas if G lies above B the equilibrium may or may not be stable. We present an alternative approach using energy principles similar to that of Erdös [E1] and Gilbert [G]. After identifying a waterline that is consistent with Archimedes' Principle, we define a unit vector normal to the waterline, by keeping the object fixed and rotating the frame of reference (waterline) by angle θ to generate all orientations satisfying Archimedes' Principle. The stable positions of a floating body occur at the minima of the potential energy. The potential energy function for a floating body is given by wheren is the unit normal vector to the waterline pointing out of the water, and θ is the rotation angle of the waterline. Note that B(θ) andn(θ) depend on θ, but G is independent of θ. Below, we derive formulas for the stable floating configurations by finding minima for U (θ).
2.2. Square Cross Section Revisited. The stability of a long floating object with square cross section, and the corresponding nontrivial floating configurations, have been investigated theoretically in a number of studies. Reid [Re] provided the first theoretical identification of stable floating equilibrium configurations based on arguments using forces and moments. Feigel & Fuzailov [FF] provided a recent alternative derivation of these equilibrium conditions, a brief review of related studies, and also detailed experiments validating the theory. Their experiments had a particular focus on floating configurations in the transition from 'flat side up' orientations to 'corner up' orientations. We revisit the square cross section configuration here with the goal of writing down the entire potential energy landscape, whose minima reveal the stable equilibrium configurations. The square has four-fold symmetry, and we exploit this in the identification of a center of buoyancy formula. A new contribution we make in the present work corresponds to situations in which the center of gravity is not at the center of the square. We specifically explore the breaking of this four-fold symmetry for floating objects with square cross sections and use the corresponding potential energy landscapes to understand the observations. In sections that follow, we demonstrate that the identification of potential energy landscapes can be obtained for shapes of more general cross sections in order to understand their stable floating configurations.
Rather than fix a waterline and consider different orientations of the square we fix a reference frame on the square with corners at (1, −1), (1, 1), (−1, 1), (−1, −1) and consider different orientations of the waterline. Three configurations are relevant as shown in Figure 2 -the first has the waterline intersecting opposite sides of the square and the second and third have the waterline intersecting adjacent sides of the square. We work out these three cases below and then give the generalization for all orientations.
For our [−1, 1] 2 square, the cross sectional area is A obj = 4. If we denote by A sub the submerged area, Archimedes' Principle requires that where R ∈ (0, 1) is the density ratio ρ obj /ρ f of the floating object to the fluid. We shall assume that the object's density is uniform throughout but if it were not the appropriate interpretation of ρ obj for the application of Archimedes' Principle would be the effective density -i.e. the object's mass divided by the volume of the object. Note that in the present context we work in terms of cross sectional area; corresponding volumes would be obtained by multiplying the cross sectional area by the length of the object in the third dimension.
Below we outline the computation of the center of buoyancy, B(θ), as a function of orientation θ for the two cases in which (1) the waterline intersects opposite sides of the square and (2) the waterline intersects adjacent sides of the square.
2.2.1. Waterline Intersects Opposite Sides of Square. Here we define the waterline by the equation where θ is the slope of the waterline and H is the y-intercept (see Figure 2). With the water assumed to occupy the region below the waterline, the submerged area can be written in terms of H as A sub = 2(1 + H). Therefore, Archimedes' Principle requires R = (1 + H)/2, or equivalently H = 2R − 1. Note that for R ∈ (0, 1) it follows that H ∈ (−1, 1). We define waterline intersection points (−1, y L ) and (1, y R ) and note that By definition, the configuration under consideration requires that y L ∈ [−1, 1] and y R ∈ [−1, 1]. Furthermore, the largest and smallest θ occur for y R = ±1, y L = ∓1, meaning that −π/4 ≤ θ ≤ π/4 and −1 ≤ tan θ ≤ 1. Combining these facts with the definitions of y L and y R , we get These can be rewritten as This range of tan θ corresponds to a range of θ values [θ min 1 , θ max 1 ] contained in [−π/4, π/4]. By symmetry, there is a corresponding configuration when rotated by ±π/2 and ±π.
The submerged area in this configuration is defined by the four points (1, −1), (1, y R ), (−1, y L ), and (−1, −1). We use (2.2) to find the center of buoyancy as the centroid of the submerged boundary region. In particular, let (x k , y k ) be given by Computing these sums, combined with the values of y L and y R and the fact that A sub = 2(1 + H), we find that the center of buoyancy takes the form where θ can take on any value defined by the inequalities (2.12). For use below we define this specific form for the center of buoyancy as B 1 (θ).
2.2.2. Waterline Intersects Adjacent Sides of Square. A similar approach can be applied to the second configuration shown in the upper right sketch of Figure 2.
Here we give expressions for the results when R ≥ 1/2 and when R < 1/2.
Case 1: R ≥ 1/2. Here we assume that the waterline intersects the left and top sides of the square at points (−1, y L ) and (x R , 1) so that three corners of the square are submerged. We consider the case R < 1/2 in the next section although note that this case can be carefully extracted from the present case (according to Gilbert [G], Feigel & Fuzailov [FF], among others).
Here we identify the waterline by This waterline cuts a triangular region of area 1 2 (1 + x R )(1 − y L ) from the original square. This means that A sub = 4 − 1 2 (1 + x R )(1 − y L ) and Rearranging this gives x R in terms of R and the waterline slope tan θ Conditions on θ come from the requirement that 0 ≤ (1 + x R ) ≤ 2 and −1 ≤ y L ≤ 1. The first of these reveals that For the case under consideration 0 < θ < π 2 . It follows that tan θ ≥ 2−2R. Equality corresponds to the waterline passing through the point (1, 1) and (−1, y L ). The other extreme corresponds to the waterline passing through the point (−1, −1) and (x R , 1). This has tan θ = 2/(1 + x R ). Here the triangular area is 1 Since θ cannot exceed this angle we have tan θ < 2/(4 − 4R). Thus, for this configuration the value of tan θ is constrained by As in the previous case, the center of buoyancy is obtained by calculating the centroid of the submerged area using (2.2). In particular, using the counterclockwise oriented vertices of the submerged polygon: we can calculate the following integral as a sum. Evaluating this integral gives For use below, we define this specific form for the center of buoyancy as B + 2 (θ). Recall that A sub = 4R and where tan θ satisfies (2.18). This condition gives a range of θ values given by [θ + min 2 , θ + max 2 ] contained in [0, π/2]. Again by symmetry, we get a corresponding set of angles by adding ±π/2, ±π.
Case 2: R < 1/2. Here assume that the waterline intersects the square at points (x L , −1) and (1, y R ) so that only the lower right corner of the square is submerged.
Here we identify the waterline by This waterline cuts a triangular region of area A sub = 1 2 (1 + y R )(1 − x L ) from the original square. This means that Rearranging this gives x L in terms of R and the waterline slope tan θ So, together these require As before, the center of buoyancy satisfies which can be calculated as a sum involving the vertices of the submerged polygonal cross section via (2.2). It follows that For use below, we define this specific form for the center of buoyancy as B − 2 (θ). Recall that A sub = 4R and where tan θ satisfies (2.27).
2.2.3. Potential Energy Expressions: Square Cross Section. As defined in (2.6) the potential energy function is given by where the unit normal to the waterline can be expressed as a function of θ aŝ n(θ) = (− sin θ, cos θ). For the square defined above with uniform density the center of gravity G = (0, 0). However, we are interested in a generalization of the square where the center of gravity, by some means, is not necessarily located at the center but rather has coordinates G = (G x , G y ). Note that as long as Archimedes' Principle is applied with the appropriate mass of the object, the calculations presented above for the center of buoyancy are independent of the location of the center of gravity. So, in what follows we treat G as nonzero in general. For ] defined by (2.12). Also define We can write the potential energy function U (θ) as follows A similar formula applies when R < 1/2 (replace B + 2 with B − 2 and the corresponding range for tan θ given in (2.27)).

Squares With
Off-Center Weights. We also explore the case of a square cross section with an off-center weight parallel to the long axis of the object. Specifically we consider 3D printed objects with a hole in the square that can be filled with a material of different density. In our experiments we had the option to leave the hole as void space or to insert a nail cut to fit the object. In either case, before floating the object tape was placed over the holes to prevent water from filling the space.
For such a configuration we can predict the modified center of gravity G = 0. In particular, consider the same square with corners at (1, −1), (1, 1), (−1, 1), (−1, −1) with a hole with circular cross section at point (a, b) with a ∈ (0, 1), b ∈ (0, 1), and radius r H . When such an object is printed there is a border around the hole whose thickness we denote by t and whose density is ρ PLA (i.e. the density of the solid print material). With the hole filled with a nail whose density is ρ nail we can compute the center of gravity of the object as a whole (printed object plus nail) as where M obj is the mass of the whole object (including the nail if one is inserted), L is the length of the object, Ω 0 denotes the cross-sectional domain of the square excluding the hole and border and Ω hole+border denotes the circular cross section that includes the (printed) border of the hole and the hole, and ρ( x) denotes the material density at position x in the plane. The square printed without a hole will have some void space in its interior and this can be controlled by changing the infill of the print. In our squares printed with a hole this infill region gets replaced by the hole plus the border material of the hole. Therefore, it is convenient to rewrite the formula (2.34) for center of gravity G as where Ω denotes the cross section of the square undisturbed by a hole. Under our assumption of a square with uniform density the first integral in this expression equates to the zero vector. That is, for a square without the off-center hole the center of gravity is located at (0, 0). This requires that the infill is sufficiently symmetric about the center of the square so that it negligibly moves the center of gravity away from (0, 0). 1 It follows that for the off-center square the center of gravity can then be estimated as where each of the density terms in these expressions are assumed to be independent of position. These integrals can be evaluated writing x = (a + r cos θ, b + r sin θ). It follows that So, for the off-center square with hole at (a, b) the center of gravity is shifted towards (a, b) from the origin by terms proportional to density differences and crosssectional areas. Note that all of the quantities in this expression can be determined by straightforward measurements and are listed in Table 1.
Practically speaking, our prints are not completely uniform in the direction orthogonal to the square face since the top and bottom square faces are solid PLA. An improved estimate for the center of gravity that accounts for the two end faces of the square of thickness t with density ρ PLA is That is, the new center of gravity, G, is shifted towards the hole location (a, b) by an amount related to the mass of the nail (we use M nail = 0 for an open hole) and terms related to the thickness of the hole and the material it replaces (either infill or boundary).
The predicted effective density for the block of square cross section is where the total mass of the object is (2.40) and V ext = s 2 L and V int = (s − 2t) 2 (L − 2t). When the block is printed with a hole of radius r H parallel to the long axis the predicted effective density is The effective density of the object without the nail filling the hole is calculated by the same formula with M nail set to zero. Various comparisons between this theory and our experimental observations and measurements are given below. First, however, we need to create our floating objects.

Methods for 3D Printing
In order to experiment with floating objects, we have opted to experiment by designing and 3D printing them. This has the advantage that we can easily create any object we can describe mathematically. In addition, we can vary the density of our print by changing the infill density, a parameter which is set at the time of printing. In this section, we describe the full workflow needed -and we have kept our presentation accessible to those with no 3D printing experience, in the hopes of making it possible for anyone with an interest to create their own experiments. The process consists of three steps: first, we need to design the objects in design software. We then need to give the print specific parameters, within a slicer software, where the choice of slicer software depends on the printer. Finally, we print the objects on a 3D printer. We detail the workflow of these steps here.
3.1. Design. We have opted to design the objects in OpenSCAD [OS], which is a free command line computer aided design (CAD) system. Figure 3 shows three of our experimental floating objects generated in OpenSCAD. Since we are interested in creating objects with a fixed cross section, we are able to do so in just a few lines of code. For convenience, we include syntax here so a reader could create their own. Copies of both sample code and stl files are available from [GIT]. To create a box with height 60 mm and a 30 mm by 30 mm square cross section, we use the command cube ([30,30,60]); For our off-center weight experiments, we have placed a hole of radius 2.5 mm lengthwise in the interior of the box. This is done by taking the set difference between the box above and a cylinder. In order to create a cylinder of height 60 mm and radius 2.5 mm, we use the command cylinder (h = 60, r=2.5, center = true, $fn=100); The command center = true centers the object so it is easier to position it with respect to the cube. The command $fn=100 indicates that the circle should be estimated by a 100-sided polygon. In order to have the cylindrical hole to be displaced from the center of the cube, we use the translate command. Putting this all together to create a vertical hole that is displaced by 10 mm diagonally from the center of the cube, we use the command sequence: ([30,30,60], center =true); translate([10,10,0]) cylinder (h = 80, r=2.5, center = true, $fn=100); } While most of our floating objects were simple cubes with and without holes, in Section 5.2 we also considered more complicated objects, where the cross section was given as a polygon in the form {(x 1 , y 1 ), (x 2 , y 2 ), . . . , (x n , y n )}. This includes for example the Mason M in Figs. 1 & 15. In order to create objects with general polygonal cross sections in OpenSCAD, we first create a polygon with the polygon command, and then we can turn this into a solid with fixed cross section using the linear extrude command. For example, the following creates a cylinder of height 60mm with a cross section which is a polygon determined by the ordered points (0, 0), (50, 0), (60, 40), (50,20). 3.2. Slicing and Printing. We now describe the printing parameters set within the slicing software. In order to speed up printing and keep objects light, a 3D printed object is usually partially hollow inside, printed with a lattice pattern filling a fixed infill fraction of the object interior. We denote this quantity by I. The fact that we can vary the infill fraction is quite useful for us, since our goal is to see how the stable floating orientations change with density. To be consistent, we did not vary the infill pattern; we printed all of our prints with the grid infill pattern. For further consistency, all of our floating objects were printed on a Makerbot Replicator 5th Generation printer and sliced in Makerbot Print proprietary software.
The mass of a print is most significantly affected by the infill fraction, but this is not the only factor. In addition, an outer border layer of the print is printed at 100% infill for a fixed thickness on the sides, and a fixed thickness on the top and bottom. For the most part, the default values for all thickness is 0.8 mm, the thickness of two shells, where a shell is 0.4 mm, i.e. the diameter of a standard extruder nozzle.
In order to calculate the predicted mass of a print, we need to know ρ PLA , the density of PLA. Most sources state that this value is around 1.24 -1.25 g cm −3 , though [DAetal] used the value 1.17 g cm −3 . This number can vary quite significantly depending on the brand of PLA, and therefore we decided to measure the value for our lab conditions. In particular, we printed a number of squares that were 30 mm by 30 mm by 1.6 mm. Since the the top and bottom of a print are 0.8 mm solid filament, the print is guaranteed to be 100% infill. We printed fifteen test tiles, roughly half with Makerbot brand filament and half with an off brand filament. Most of our test tiles were printed in part on the same Makerbot printer and slicer as the floating objects, but we also printed some on a Monoprice Mini. We found that the printer used was not an important factor for the mass of the test tiles, but the brand was an important factor. For Makerbot brand filament and off-brand filament respectively, our measured densities were 1.082 ± 0.017 and 1.152 ± 0.010 in g cm −3 . We have chosen the value ρ P LA = 1.15 g cm −3 for our calculations.
Taking into account the infill percentage, shells, and measured value for ρ P LA , we expect the mass of our object to be given by ρ PLA times V ext , the volume of the outer layer plus ρ PLA times I, the infill fraction, times V int , the volume of the interior of the print. If there is a nail embedded in the hole, then we additionally have to include this value in our calculation. We have described this calculation above in (2.41). Figure 4 shows the predicted values for prints without a nail hole (dashed red) with a nail hole (solid green). We also show the measured values for the prints we have made without a nail hole (red x) and with a nail hole (green o). Note that while quite useful to be able to predict mass of a print from the infill density, all of our calculations of R that we used for particular prints in our experiments are based on the mass of the print measured directly using a digital scale.
It is critical for the print to be uniform in the longest direction, as our calculations consider this direction as completely invariant. Therefore we printed with the long direction going from top to bottom so that the infill is identical at each cross section.  Since our prints were printed on a Makerbot printer, we have used the Makerbot proprietary slicing software Makerbot Print. However, there is nothing about the process that could not be modified to other slicer software.
The objects are ready to float!

Experiments and Data Acquisition
Once the objects were printed our objectives turned to floating these objects and obtaining measurements associated with the stable floating orientations that we could compare with the theory. After floating the objects in a test tank to observe qualitative behavior (e.g. confirming the objects floated, identifying the number of stable orientations, etc.) we turned to gathering quantitative information in the form of stable floating orientations.
We used a tank with clear flat sides and filled with tap water (see Figure 5). In order to measure angles associated with a stable floating orientation the object was carefully placed by itself in the tank and allowed to come to rest. In order to ease the task of keeping the object in a position with its square cross section facing the side of the tank where the camera was positioned we placed vertical guides a small distance away from the object. One of these guides was a ruler that we could later use in the image analysis for calibration. A typical view of a floating object is shown in Figure 6.
Digital images were obtained with a Canon EOS Rebel XSi used in manual focus mode. A remote shutter release was used to avoid bumping the camera, which was mounted on a tripod sitting on the same table as the tank. Care was also taken to line up the camera at the water level and to have the object floating Figure 5. Our experimental setup involves floating 3D printed objects mostly with square cross sections. An occasional print had density ratio R > 1 and did not float. Figure 6. A typical view of the experimental set up to measure angles associated with stable floating orientations. Vertical guides were used as an aid to keep the object floating in the part of the tank on which the camera was focused. Care was taken to assure that contact with underwater objects was avoided and any incidental contact of the object with the guides did not alter its floating orientation. We also took care to prevent contact of object with the sides of the tank.
in a direct line to the camera. A certain amount of refraction could be observed, especially in directions off the main viewing axis.
We used Matlab's grabit.m software to extract information from each image such as the one in Figure 6. We first calibrated the view by selecting points on the ruler in view on the right, taking care to use points above the waterline to avoid distortion of the ruler by the water. For a floating square we then identified 6 points in the image. Two of these points were chosen on the waterline relatively far from the object on both sides so that the orientation of the waterline could be obtained. Four other points were chosen at the four corners of the square (clockwise starting from the top of the square). With this information we obtained four vectors by taking differences of adjacent vertices and used cos θ = a T b/ a 2 b 2 where b is a waterline vector and a is one of square-side vectors. In the example of Figure 6 one pair of sides (northwest and southeast) were used to identify an estimate for θ while the other pair of sides (northeast and southwest) were used to estimate the complementary angle with respect to π/2. In the case of the Mason M, we used a similar approach but chose four points along the bottom two 'legs' of the M to make angle measurements.
For prints with holes, we covered the holes (void or filled with a nail) with waterproof tape to keep water out. In cases with no hole or when the hole was in the center of the object, we used a very small amount of nail polish to mark one corner of the square. This allowed us to assess asymmetry of the object that was present either unplanned or by design.
For each object we first measured its mass (with hole either left as a void or filled with a nail) using a digital scale and from that obtained an effective density according to ρ obj = M obj /V obj . This corresponds to the material density for a uniform object or the effective uniform density for a non-uniform object. The volume of the object V obj was computed either from measurements of the dimension and known formulas (e.g. width times length times height for a rectangular box) or using polyarea in Matlab (for more general cross sections) multiplied by the 3D print scaling.

Computational Approaches
Various Matlab codes were developed to compute results and analyze our floating shapes.
5.1. Square Cross Sections. Various Matlab codes developed for the square cross sections have been posted in a GitHub repository [GIT]. These include • SQUARE PE GxGy.m: This code is based on the potential energy formulas outlined in the section on the square. It generates for given values for the density ratio R and the center of gravity (G x , G y ) the computed potential energy landscape and includes options to plot the potential energy landscape and the square floating in a stable orientation. This code was used to generate the theoretical predictions in Figures 8, 9, and 10, for example. • SQUARE ANGLES GxGy R Looper.m: This code is based on the formulas given in the section on the square and plots for a given center of gravity (G x , G y ) and specified range R ∈ [R min , R max ] the stable orientation angles. This code, for example, was used to generate Figures 7 and 13.

General Polygonal Cross Sections.
In the case of a long floating object of uniform density with a general polygonal cross section, though it is no longer possible to give as detailed an analysis as in the case of a square, we are still able to apply Archimedes' Principle and calculate the center of gravity, center of buoyancy, and the potential energy of a floating configuration, as we describe below. 5.2.1. Computation of Stable Floating Configurations. When doing the calculations for a general polygon, we wrote a program that takes two vectors with the x and y values of our polygon and a density ratio and goes through the following algorithm. Firstly, it takes the shape of the given polygon and calculates the center of gravity of the object, assuming a uniform density throughout. Secondly, we identify the correct placement of the waterline determined by Archimedes' Principle that establishes the correct submerged area to total area ratio. This is done by use of a bisection method where at an orientation we take the lowest point of our object and create a waterline through it and make that our lower bound. We then take the highest point of our object and make that our upper bound. We then calculate the area ratio of each of our bounds and find the midpoint of our upper and lower bounds and calculate the area ratio for the waterline going through the midpoint. If the waterline through the midpoint is above the correct placement, it becomes our new upper bound and if it is below the correct placement, then it becomes the new lower bound. We apply the bisection method until the correct waterline is found for our original orientation. The correct placement of our waterline in the case of uniform density is where the waterline splits the object into two areas where the submerged area relative to the total area is equal to the desired density ratio (e.g. for an iceberg with uniform density, the line will be such that the new submerged area created by the line relative to the area of the original polygon will match the density ratio 0.8912). Thirdly, we compute the center of buoyancy of the object by applying the same method used for the center of gravity except using the submerged area determined by the polygon and the waterline. Finally, we calculate the potential energy function for our polygon at the orientation and repeat the process for all angles.
The input to this code is a planar polygonal region, oriented counterclockwise and a density ratio of the object relative to the water. The output is a plot of the potential energy landscape with respect to the angle. By default, the computation is performed for uniform density. However, it is possible to compute this information for objects with non-uniform density if one inputs the center of gravity.
6. Results 6.1. Floating Squares: Symmetric Case. Figure 7 shows stable equilibrium angles as a function of density ratio R for objects with square cross sections and center of gravity at the center of the square, G = (0, 0). Various experimental results are shown for 3D printed shapes with different effective densities. These effective densities have been modified as described earlier by adjusting the infill as well as printing objects with a hole at the center which we can leave as void space or fill with a denser object, such as a nail.
As a visualization, we show several potential energy landscapes and selected shape orientations predicted from the theory and observed experimentally. Figure 8 shows the potential energy landscape for G = (0, 0) for a case with R = 0.23296 which corresponds to a region in parameter space where eight stable orientations exist. The eight orientations in this case come in pairs, as indicated in the lower portions of Figure 8. The eight experimentally-observed orientations are shown by the points at R = 0.23296 in Figure 7. 2 Figure 9 shows the potential energy landscape for G = (0, 0) for a case with R = 0.4856 which corresponds to a region in parameter space where four stable orientations exist. These orientations correspond to the object floating with the corner straight up. Experimental measurements for these angles correspond the angle measurements shown at R = 0.4856 in Figure 7. Figure 10 shows the potential energy landscape for G = (0, 0) for a case with R = 0.9322 which corresponds to a region in parameter space where four stable orientations exist. These orientations correspond to the object floating with the flat side of the square straight up. Experimental measurements for these angles correspond the angle measurements shown at R = 0.9322 in Figure 7.
A more thorough experimental exploration of the parameter space, particularly in the region around R = 0.25, where eight stable orientations can be identified has recently been done by Feigel & Fuzailov [FF]. Those authors used a larger floating object (114 mm × 114 mm × 353 mm) constructed using two tin tea boxes that could be fitted/weighted with additional bars and magnets to adjust the object's effective density.
Figure 7. This plot shows stable floating orientations versus density ratio R for G = (0, 0). The marks show various measured equilibrium orientations for several of our 3D printed objects.
6.2. Floating Squares: Breaking Symmetry. Figure 11 show the potential energy landscape for three different prints corresponding to square cross sections with nail-filled holes at A:   Figure 9. The upper left plot shows the potential energy landscape for the square with G = (0, 0) and R = 0.4856. There are four stable equilibria corresponding to θ = π 4 ± π 2 n for integer n. The upper right plot shows that the square floats with vertex pointing upwards. These four stable orientations are also observed experimentally (one such orientation is shown in the image). Measured angles for this case are shown in Figure 7.
in units of s/2 where s is the length of the side of the square. These three have density ratios of R = 0.4856 for Case A, R = 0.4874 for Case B, and R = 0.4911 for Case C.
For Case A with nail-filled hole at (0, 0) already discussed in Figures 7 and 9, four stable orientations, with a corner of the square pointing straight up, are predicted theoretically and observed experimentally.
Experiments for Case B, with a nail-filled hole at (0.3, 0.3), are shown in Figure 12. In this case the center of gravity G is no longer at the center of the square and the symmetry is broken. Despite this broken symmetry, four stable orientations are still observed experimentally as shown in Figure 12. Our theory predicts that (G x , G y ) = (0.06789, 0.06789) and the corresponding potential energy plot is shown by the solid curve in the upper right plot of Figure 11. There are only two stable orientations predicted at this value of G and so our theory does not match the experimental observations. However, two other curves are shown in the upper right plot in Figure 11 -the dashed curve has (G x , G y ) = (0.05, 0.05) and the dash-dotted curve has (G x , G y ) = (0.04, 0.04). These curves with nearby values of  Figure 10. The upper left plot shows the potential energy landscape for the square with G = (0, 0) and R = 0.9322. There are four stable equilibria corresponding to θ = 0 ± π 2 n for integer n. The upper right plot shows that the square floats deep in the water with flat side up. These four stable orientations are also observed experimentally (one such orientation is shown in the image). Measured angles for this case are shown in Figure 7.
G indicate that there are two other stable orientations nearby. One explanation for the discrepancy between experiment and theory is that the center of gravity of our print is not exactly where we predict it to be, perhaps due to uncertainties in the infill structure of the print. Another potential source of imprecision in the center of gravity is the nail not fitting precisely in the middle of the hole, but the nail seems to fit snugly in the hole, so it seems a less likely explanation. Also, there is evidence of menisci at the solid-liquid-air contact line in Figure 12 suggesting that surface tension could provide a large enough force to hold the print in an otherwise slightly unstable configuration.
For Case C, with a nail-filled hole at (0.45, 0.45) and (G x , G y ) = (0.1017, 0.1017) the potential energy landscape shown in the lower plot of Figure 11 indicates that only two stable equilibria exist. For this case both theory and experiment are in agreement on the number of stable equilibria. The experimental images for this case are not shown but are very similar to the upper left and lower left images of Figure 12. Corresponding stable orientations with the nail to the left or right as in the upper right and lower right images of Figure 12 no longer exist for Case C. The open red circles indicate the theoretical local minima of the potential energy. The vertical red dashed lines show the angles at which experimentally-floating squares appear to be stable. Case A has four-fold symmetry and four stable orientations are predicted theoretically and observed experimentally (see also Figure 9). For Case B, in which the symmetry is broken, our theory predicts only two stable orientations but we observe four experimentally (these four orientations are shown in Figure 12). The dashed and dash-dotted black curves in the upper right plot show two other potential energy landscapes for nearby values of G (see text for details) indicating the presence of nearby stable states. For Case C, the square is farther from symmetric and both theory and experiment predict only two stable orientations. We do not show experimental images for Case C, but they are similar to the upper left and lower left images in Figure 12.
As we have just observed, when G = 0 the symmetry of the square is broken and the predicted number of stable equilibria change. We display this in Figure 13 for four cases G = (0.01, 0.01), (0.02, 0.02), (0.04, 0.04), and (0.06789, 0.06789). On the last two of these we have plotted the four experimentally-observed stable orientation angles for the print denoted Case B with a nail-filled hole at (0.3, 0.3).  Figure 11 with a nail-filled hole at (0.3, 0.3) giving a nonzero center of gravity G. When the hole is up or down the stable orientations correspond to corner straight up. When the hole is left or right the top corner of the square is deflected counterclockwise or clockwise, respectively.
As just discussed for Case B, our theory predicts G = (0.06789, 0.06789) and only two equilibrium values of θ exist for this case as shown in the lower right plot of Figure 13. Reduction in the value of G to (0.04, 0.04), for example as indicated in the lower left plot of Figure 13 shows four stable orientations at the density ratio, R = 0.4874, of Case B. The upper two plots corresponding to G = (0.01, 0.01) and (0.02, 0.02). Especially in comparison to the predicted stable angles shown in Figure 7 these plots reveal interesting bifurcation structure as the square symmetry is broken. Note that in addition to floating squares with four or eight stable orientations, when the symmetry is broken in this way (moving the center of gravity towards a corner) there are also situations where either three or six stable orientations are predicted. pixel to millimeter (mm) conversion factor used to scale for 3D printing gave the area. In particular, we found that A M = 126, 550 · (0.06) 2 mm 2 = 455.58 mm 2 . The length of our Mason M was 70 mm, giving a volume of 31, 891 mm 3 . Since its mass was 27.92g its corresponding effective density is ρ eff = 0.8755 g cm −3 , which, in comparison to the density of water, is fairly close to the typical density ratio of an iceberg in a polar sea. Figure 14 shows the potential energy landscape along with predicted stable configurations for the Mason M. The four stable orientations of the floating Mason M are shown in Figure 15. Figure 16 shows another variation of the full Mason M (with feathers) inspired by calving events that can happen with real icebergs. These two plots show stable floating orientations, computed using the general polygonal code, where a calving event of shedding a feather results in a new stable floating orientation (right plot). The red line in the plot on the right shows the pre-calving waterline position. Results on their stability associated with metacentric concepts that date back to eighteenth-century continue to play a major role in areas of naval architecture and the study of icebergs.
The field has recently gained mathematical interest, especially on Ulam's floating body problem. In problem 19 from the Scottish Book Ulam asks: "Is a solid of uniform density which will float in water in every position a sphere?" [M]. Counterexamples to Ulam's problem in the plane for ρ = 1/2 go back to Auerbach [Au], and for ρ = 1/2 to Wegner who also obtained results for non-convex bodies (holes in the body are allowed) in [W2]. Most recently Florentin et al. [F] gave an affirmative answer for a class of origin symmetric n-dimensional convex bodies  with ρ = 1/2 relative to water. While counterexamples in the plane exist, Ulam's floating body problem in higher dimensions is open to the best of our knowledge.
7.2. Limitations. Using 3D printed models we were able to produce and collect experimental results in agreement with the theory. Our results and experimentation process can be reproduced readily, increasing access to 3D printers within educational facilities gives way for students to experiment with stability and generate different scenarios including Ulam's floating body problem.
The biggest challenge faced when producing our 3D samples is related to density. To design infill patterns and infill densities using computer-aided design(CAD), we produced samples to measure precision and increase the accuracy of our experiment and identify possible sources of experimental errors. A limitation in the 3D printing process is the difficulty in producing low-density objects to desired accuracy. To maintain an even distribution of mass, considerations must be taken with respect to the infill pattern. Zigzag, grid, triangle, and concentric patterns are recommended when an even distribution is desired. While Makerbot's cat infill model is not recommended for experimenting with homogenous objects as it is asymmetric, it made for interesting observations and ultimately motivated further investigation into non-homogeneous objects, and it's also really cute, cf. Figure 17.
Our work here has focused on effectively two-dimensional shapes. Mathematical challenges for 3D floating shapes have been examined (e.g. Erdös et al. [E2] and Wegner [W1]) and of course most realistic floating objects such as icebergs are three dimensional. Approaches using 3D Print design are likely to prove highly useful for future studies in these directions.