Examples of Non-Uniqueness of the Equilibrium States for a Floating Ball

We provide a numerical algorithm for numerically approximating a centrally located floating ball. We give examples of equilibria, and we present non-unique cases for the same physical parameters when the density of the ball is either greater than the supporting liquid (heavy) or lighter than the density of the vapor above (light). We classify the non-uniqueness by analyzing a function related to the force balance. We derive the potential energy of these states, and make comparisons of the non-unique cases. In the cases of both the light and heavy floating balls, the evidence presented supports the conjecture that when there are two equilibria, the one with lower energy corresponds to the location of triple junction (between the ball, the vapor and the liquid) that is closer to the equator of the ball.


Introduction
Consider a ball of density B ρ floating at the surface of a fluid that has density ρ  .We present numerical examples of non-uniqueness of the equilibrium states for these configurations.We will also provide a framework for the classification of these states, including an energy analysis.The energy analysis is used to determine which of the two equilibria has the lower energy, and thus, at least amongst centrally located floating balls, this process finds the energy minimizing configuration.Under these conditions, this is the configuration that our model predicts which will be found in experiments.We begin with a precise formulation of our model.
The energies considered in this model are due to the surface tension and gravity.Surface tension energy is taken, as usual, to be proportional to the area of the free surface with proportionality constant σ called the surface tension constant.The energy due to gravity consists of two terms, one corresponding to the liquid and another to the floating object.The former is proportional to the density of the liquid, ρ  , a gravitational constant, and a volume integral of the physical height, z, in the gravity field.The latter is similar except the density of the floating object, B ρ , is used, and the integral is taken over the volume of the floating object.If we denote the floating ball by B, the liquid by E, and Λ the free liquid-air interface, w  the wetted portion of the bounding walls, if they exist, B  the wetted portion of the floating ball.Then the energy of this configuration is given by d d , with wetting coefficients = cos γ and w γ at the contact with the ball and the wall, respectively.Here z is height in the vertical direction, and dV is the volume measure.Also, ⋅ denotes surface area, calculated using a Hausdorff measure, as appropriate.It is often con- venient to refer to the mathematical energy of the system, which is merely the scaled energy σ  .Note that then the gravitational terms become d and d with "capillary constant" g κ ρ σ =  .The natural physical setting for these problems is in a bounded container in 3   .We will study the lower dimensional problem here for two reasons: first, there is a growing literature of flotation in this setting, and second, there is an approach to prove the existence (and, when applicable, uniqueness) of the equilibria in this setting using a phase plane analysis; and it is useful to have robust numerical simulations.We will consider the settings in 2   of both a bounded container and an unbounded sea of fluid.We can interpret these problems as the lower dimensional setting, or we can interpret them as unbounded in one horizontal direction.In the latter, the floating ball may be seen as an infinitely long log floating in an infinitely long trough.The triple contact line of the liquid with the air at the surface of the floating object is then a straight line extending along the unbounded log.In what follows we preserve the intuitive meaning of area, volume, and energy by interpreting the lower dimensional setting as generating cylindrical configurations in the unbounded horizontal direction, and then we take the area, volume, and energy to be per unit distance in that unbounded direction.
The goal of this note is to provide numerically computed approximations to some sample configurations when the ball is centrally located in the container and the fluid is assumed to be symmetric about the central axis.We will focus on the cases of the density B ρ outside of the range ( ) 0, ρ  .One could classify the admissible densities into three types: light, medium, and heavy.This would correspond to ( ) , respectively.Specifically, we will classify the behavior of a certain function

( )
F φ that is pi- votal in the force balance considerations, as derived using variational techniques.Then we will follow this by a measurement of the mathematical energy, specifically of interest for the cases of non-uniqueness.For a theoretical treatment of this setting in general, see first McCuan and Treinen [1], and also McCuan [2] for further details.We will use the results of those variational arguments in what follows.
With this setting established, we are able to proceed to the configurations of interest.Assume that the ball floats centrally, and so that the center of the ball is at a height d.Thus the boundary of the ball can be described by its azimuthal angle φ and its radius a, given a height d.The fluid-air interface contacts the ball at a particu- lar azimuthal angle, denoted by φ .In the bounded cases a volume of fluid is described by the wetted container walls, the fluid-air interface, and the wetted surface of the ball.This volume is held fixed, and introduces a Lagrange multiplier λ into the problem.The height of the fluid is given by the differential equation where 2H is the mean curvature operator, commonly seen as when the surface is a graph over some base domain.However, we will not restrict ourselves to these limited configurations.We have boundary conditions cos Tu where ν is the exterior unit normal on the container, or the interior unit normal on the ball, as appropriate.Al- so, γ is taken to be the contact angle along the container wall, or the ball itself, also we assume that γ is con- stant, up to possibly having different values on the container wall and on the ball, as described above.This restriction that γ is locally a constant is simply a statement that we are considering only uniform materials in this current work.The underlying model is flexible enough to accommodate non-constant γ , if γ is well behaved.
In [1], we derived an additional necessary condition that does not appear in the classical literature consisting of fluid interactions with rigid solid objects.A manuscript by Finn [3] is the standard reference for the classical literature.We need to define several objects before we can state this condition.Denote the volume of fluid by  .Set n  to be the outward pointing unit co-normal along the boundary of the liquid-air interface.Set N to be the unit normal to ∂ pointing out of  .The fluid-air interface has a surface tension σ , and the potential energy due to gravity is measured as where G depends on position and the material.For our application, if z measures height in the vertical direction, then G gz ρ

=
, for the appropriate density ρ .Then the con- dition for free floating is Under the assumption of symmetry about the vertical axis, the PDE with boundary conditions can be reduced to This representation of the equation allows for the parameterized curve to pass through both inflection points and vertical points.Various authors have studied this system with differing boundary conditions and also as a family of solution curves.Solutions are sometimes known as Euler elastica.See Aspley, He, and McCuan [4], Euler [5], Giaquinta and Hildebrandt [6] (pp.142-144), McCuan [7] and [8], and Wente [9] for both historical origins, as well as applications to capillarity.
Then (2) becomes and we define ( ) F φ to be the left hand side of this equation.A key observation is that the locally observed information is collected on the left side of this equation, and the right side contains the global density information.In [1] we interpreted this as a version of a force balance condition.We seek solutions to (3) that satisfy (4) and where 0 0 v > is some prescribed quantity of volume large enough that the ball need not touch the container.This problem is the bounded container problem in 2   .The unknown quantities are , , , d φ λ  .The details are in Section 2.1.If we replace the bounded container with an infinite sea of liquid, then we follow Bhatnagar and Finn [10] in taking the Lagrange multiplier λ to be 0, and the condition π 2 B ψ γ = − when ( ) . We then seek solutions to (3) thus modified that satisfy (4) with 0 λ = .This problem is the unbounded container problem in 2  .This problem is significantly simpler, and one may view it as adding assumptions compared to the bounded container problem, the result of which is a family of curves representing the fluid interface that is completely characterized.In fact, it can be shown that not only is there existence and uniqueness of the boundary value problem in this setting, but we also have a formula for both ( ) x ψ and ( ) u ψ .(See [10], though unfortunately the formula is stated incorrectly there.)Thus the unknown quantities in the 2D unbounded container problem are reduced to φ .The details are in Section 2.2.

Methods
We will consider first the bounded container in 2   , followed by the unbounded container.In both of these con- figurations there is a common method of approach.The floating ball problem is cast as solving a system of ODEs coupled with side conditions for the volume constraint and the new free boundary condition.The side conditions can be considered simultaneously with the standard boundary conditions, and form a set of necessary conditions.We employ shooting methods for the underlying systems of ODEs coupled with nonlinear zero finding algorithms for vectorized necessary conditions.The zero finding algorithms require initial guesses for the free parameters, which we tune to the given physical with estimates when we are able.
Next we illustrate how the new free boundary condition can be used to determine when there is non uni-queness of the equilibria.We give numerical examples.Finally, in the case of the bounded container, we calculate the potential energy of each configuration and compare the energy of the two configurations with the same physical properties.

Bounded Container in  2
We first need to measure the volume of fluid held in the container.One may compute ( ) ( ) ( ) Our formula does not require that the height u is a function over a base domain.
In what follows we will assign 1 κ = in all of the numerical computations.This can be seen as applying a standard scaling argument, however we are then unable to further normalize the container described by R, and the radius of the ball a.Thus these two parameters are inherent features of this model.
The numerical procedure is a shooting method, where we integrate the system of ODEs cos sin to some ending arc-length  using ODE45 in Matlab.We use the finest permitted tolerances, with the absolute tolerance set to 1 14 e − and the relative tolerance set to 2.23 14 e − .We have free parameters , , d φ λ , and  .
We use these to satisfy the physical conditions which say that the configuration satisfies the volume constraint and the force balance condition as well as the need for the fluid interface to extend to the wall of the container with the prescribed contact angle.The parameters , , , d φ λ  are given the following initial guesses: This leads to a solution of the initial value problem where we can evaluate the conditions ( 8)- (11).We use Matlab's fsolve to then vary , , , d φ λ  , computing the solution to the IVP out to the arc-length  at each step, until Equations ( 8)- (11) are satisfied up to the requested tolerances.Matlab's fsolve uses a dogleg variant of a trust region method to solve this 4-dimensional zero finding problem.Here we use tolerances that are coarser than the underlying ode solver's tolerances, as the zero finding depends on those approximations.Specifically, we use termination tolerance on the objective function value of 1 12 e − and on the input values of 1 12 e − .We next proceed with a few examples.In Figure 1 we set 2 R = ,  .This leads to a solution of the IVP which we can use to evaluate the following conditions to within a given tolerance: ) ( ) We then calculate the gravitational potential energy of the ball: The gravitational potential energy of the liquid is more convenient to compute in sections.First we will treat the portion directly below the interface, but due to the lack of a closed form expression and also due to the variable step size of the data we do this numerically with the trapezoid rule:    where 0 sin x a φ = .The error here is bounded by where h is the maximum step size taken in the irregularly spaced data, and ξ is some number in the interval , then we simply need to add to this quantity the portion directly under the ball that has not yet been accounted for: ( ) The case where π 2 φ < has the complication that we must remove the portion of the fluid where the ball overlaps the portion directly below the free interface.We find the formulas match if we additionally compute the energy from the top of the ball, then subtract the entirety of the corresponding energy of the ball: which matches (15) upon subtracting 2 πa d , the energy of the whole ball.

Unbounded Container in  2
We next consider the unbounded container problem in 2  .This is envisioned as an infinitely long log floating on an unbounded sea of liquid.As such, there is only the contact angle on the ball, so we use B γ γ = .It is worth noting that the only relative relation remaining is that comparing a to κ .In the case of the bounded container there is the significantly more complicated interaction between a, R, and κ , as well as the contact angle w γ .
We are first interested in finding a solution to cos sin = .This can be explicitly integrated: ( ) ( ) ( ) See [10] and [9].We proceed to the full problem of computing the floating ball configuration in the unbounded 2  configu- ration.Using our formulas, for each to fix the height of the ball.We are then able to evaluate the necessary condition with 0 λ = .We are left with this single equation as well as one free parameter: φ .We use Matlab's zero find- ing function FZERO to vary φ , computing the solution to the asymptotic boundary value problem at each step, until this equation is satisfied to within given tolerance.
Our results can be compared to Bhatnagar and Finn [10], where an independent variational formulation was used.In our case, we are formally applying the results from [1] which assumed a bounded container in order to measure the energy of the configuration, whereas in [10] their variational argument was constructed with care to account for the difficulties of infinitesimal variations of infinite quantities.

Results and Discussion
We begin here with the floating ball in a bounded container, and we ask the It should be stated that while this current work focuses more directly on the influence of the contact angles on the nonuniqueness criterion, the influence of R and a on the nonuniqueness criterion is just as important.To see this impact, compare Figure 9 with Figure 10.We leave a separate study focusing on fixed contact angles and varied values of R and a for a further work.
With the energy computed as in Section 2.1, we first use it in conjunction with the graph of ( ) we compute the mathematical energy of that configuration.This curve is plotted with that of ( ) F φ in Figure 11.Carefully note that the density of the ball changes along this curve.
We do not specify a fixed density, however, we back out the density that appears in the graph of ( ) Then the energy is computed with this value of B ρ .Further, the energies of the two equilibria are compared in Figure 11, and at least for the particular height picked to illustrate this energy difference, the equilibria with   and for each of these, we use the value of ( ) F φ to interpolate the smaller value of φ that shares the same height of F. Finally, we interpolate the energy at this smaller value of φ .
Then we plot the energy for φ in ( ) and the corresponding energy values that come from this comparison scheme.The results are shown in Figure 12.The corresponding comparison is done for the heavy floating ball, see Figure 13.In both of these cases the value of φ further away from the endpoints of the interval ( ) 0, π had lower energy than the other possibility.This is some evidence for the conjecture that the energy mi- nimizer contacts the ball closer to the equator when compared to other equilibria.We will state this conjecture more precisely below.
We then included the above analysis into the program that generated Figure 9 with a 50 × 50 grid of values for B γ and w γ , and the results are shown in Figure 14.A series of these simulations were performed amounting to 3.5 million tests of the conjecture.The results are in Table 1.The cases with 3 R = and 1 a < were attempted, however the initial guess for the zero finding function needs to be refined in these cases, as those guesses had been tuned to configurations where the curvature was larger, and convergence to semi-equilibrium is not uniform on the grid considered.The cases 2 R = , 1.75 a = and 3 R = with 2 a > are at the limits of the ability to prescribe finer tolerances for the problem and the results are not conclusive for the cases we have considered there.We should note that in the case 2 R = , 1 a = we needed to prescribe the relative tolerance of Matlab's ode45 to 2. 23 14 e − in order to verify the conjecture there.This became the standard for our tests, and as the spacing of floating point numbers is 2.2204 16 e − , we are currently unable to achieve finer results.
Next, we trace the energy values from a sampling of balls with given, fixed densities as the semi-equilibrium are computed through the range ( ) . This is in contrast to the energy graphs that we have so far   Table 1.The configurations where the conjecture was tested, with the displayed radius a and the half-container width R.
Each entry on this table represents 250,000 configurations with 50 × 50 evenly spaced points of ( ) examined.Here we fix a few densities and in contrast, in the previous discussions the densities would vary with the value of ( ) 15 shows the results of this, which should be seen as following a ball as it rises from the depths, first registering when it contacts the interface with 0 φ = in a semi-equilibrium state, and passing continuously through that parameter until leaving the fluid at π φ = .It should be understood that this is not a continuous behavior in terms of the center height d.In fact, in order to achieve the semi-equilibria, in general d changes discontinuously at both the point of beginning and ending the contact with the free interface.We do not include that behavior in our graphs.The feature exhibited in these figures that that there are at most two interior local extrema along the trajectory of the energy profile curves.This is in agreement with the maximum of two equilibria for a fixed density, of either a light or heavy ball.The phenomenon can be described at follows, as the ball passes through a full range of heights, though not necessarily monotonically.First, let 0 B ρ < .Then, as d increases from the depths, the energy decreases.At the contact with the interface there will be a discontinuous change in d.Then we move through the parameter φ , increasing from 0 φ = as the energy continues to decrease.Then there appears a local energy minimizer, where the surface energies play a stronger role than the gravitational energies.From this point the energy may increase from a local minimum, or decrease from a saddle point.On the remainder of the curve there may be a second equilibia, which is a local energy maximum.The energy is continuous to the endpoint where π φ = .At this point there is another discontinuous change in the height d as the ball frees itself from the interface, and from this point it continues upwards indefinitely as the energy continues to decrease.Second, let B ρ ρ >  .Here the ball rises from the depths, with increasing energy up to the point of contact with the interface.Then there is a discontinuous change in d at the contact, but where the configuration is in a semi-equilibrium at 0 φ = .Then as φ increases, the energy decreases to a local energy minimizer.From there the energy increases, and there may be another equilibrium along the curve, which would be a local energy maximum.Either way, the energy is continuous up to the point where π φ = , and then there is another discontinuous change as the ball leaves the liquid.From there the ball continues upward, and the energy increases with d.
Finally, we turn to a study of ( ) F φ for the floating ball problem in an unbounded container.Adapting the argument from the case of an bounded fluid in 2  , we generate F over [ ] 0, π .In this setting, however, once φ is removed from the list of unknowns, and  is removed from the necessary conditions, we are able to generate a configuration for each φ using the explicit solutions as above. .The behavior is in line with the simulations done with the bounded case in 2   .We do not here proceed with an energy comparison between these non-unique cases for two reasons.First, the one approach that might be used would be to fix some large domain with which to measure the energies and compare the two values obtained.This is seemingly somewhat arbitrary, and has the possibility of not detecting the correct case when the comparison is close.Second, and more important, is that Bhagnagar and Finn [10] carefully analyzed the energy comparisons, and their method is superior to that just described.We do not wish to repeat their analysis, so we simply take note that their analysis supports the following conjecture in the examples that they computed explicitly.Further, Bhatnagar and Finn found examples of nonuniqueness, and their energy based methods give rise to a more   rigorous stability analysis.The results from this section then can be seen as purely an analysis of the non-uniqueness criteria.Conjecture 1.The centrally located floating ball has at most two equilibria in a bounded container in 2D, and the centrally located floating ball has at most two equilibria in the unbounded configuration in 2D.In both cases, this will occur at some value or values of the azimuthal angle φ on the ball.In the case that 0 B ρ < , if there are two configurations then the solution with the smaller value of φ has lower energy.In the case that B ρ ρ >  , if there are two configurations then the solution with the larger value of φ has lower energy.

Conclusion
We have developed a robust numerical solver for finding the equilibria of a centrally located floating ball in both bounded and unbounded problems in 2  .The non-unique cases for both the light ball and the heavy ball have been analyzed using the function ( ) F φ , and the uniqueness and non-uniqueness depend on the geometry of the graph of that function.We calculated the potential energy of the bounded configurations, and used that first, to determine which non-unique equilibrium was the local energy minimizer, and second, to trace the energy profile of a ball as it is passed through the fluid interface.Our methods were used to formulate a conjecture on the energy minimizer, and 3.5 million test cases were verified.

γFigure 1 .
Figure 1.Comparing contact angles w γ set to π on the left and 0 on the right.

Figure 2 .
Figure 2. Non-uniqueness: two values of φ for the same physical parameters.

Figure 3 .
Figure 3. Computing partial solutions for values of φ .

Figure 4 .
Figure 4. Computing partial solutions for values of φ .

Figure 5 .
Figure 5. Symmetric values of contact angles.
and lim s x →∞ = ∞ .The curve ( ) , x u forms a generating curve for the cylindrical liquid-air interface.This system is equivalent to the following system

.
The configuration on the left is conjectured to be stable, with a value of 0.5732 φ =, and the configuration on the right is conjectured to be unstable, with a value of 0.14686 φ =.In comparison, there was only one solution when ( )

φφ
following question: what values of F as above and then test the regions 2π F > and 0 F < for monoto- nicity.The results appear in Figure 9 with Figure 10.It is worth mentioning that there are grid points on each axis, and( ) F φ was evaluated at 100 grid points of φ , giving 250,000 numerical solutions to the (par- tial) floating ball problem for each of these figures.The cases marked with a star represent curves that begin at 0 to its global maximum, from which it decreases to its terminal value at π only.The cases marked with a plus represent curves that begin at 0 - ceeding to its global minimum, from which it increases to its terminal value at π

Figure 11 .
Figure 11.Using the function F to determine which solution has the lower energy.The curves are F and the energy.The lower horizontal line fixes a height that picks up two values of φ for the value of F (and thus B ρ ).These two φ values are indicated by the vertical lines.Finally, the upper horizontal line is shown that meets the intersection of the energy graph and the rightmost vertical line.Note that the energy is below the intersection of the upper horizontal line and the left vertical line.This implies the configuration with the smaller value of φ has lower energy.

Figure 12 .
Figure 12.The light ball energy comparison.

Figure 13 .
Figure 13.The heavy ball energy comparison.

Figure 15 .
Figure 15.Given 5 sample densities B ρ , tracing the energy of the system as φ moves from 0 to π .On the left

Figure 16 . Figure 17 showsF
shows the results for 0 γ = and π γ = .Again, note the lack of a possibility of a solution existing if B an example when it is possible to have non-unique solutions for both of the cases B < for monotonicity.The results are collected in Figure18, where we display the results for a selection of radii Figure 16.( ) F φ with contact angles 0,π γ = .

Figure 17 .
Figure 17.An example of ( ) F φ with nonuniqueness for both light and heavy floating balls.