# Iceberg dynamics

Prompted by this tweet and campaign for icebergs to be depicted in their most stable equilibrium orientation, here is a Python script modelling the dynamics of a two-dimensional iceberg which starts in an arbitrary orientation and position and relaxes under gravitational and buoyant forces to its most stable configuration. A cork floats "on its side": with its longest axis parallel to the water's surface (it doesn't bob around with its longest axis vertical), and an iceberg does the same.

The centre of mass ($G$) and the centre of buoyancy ($B$, the centre of mass of the displaced water) are indicated by the red and blue dots respectively. The motion of the iceberg (treated as a rigid body) is dictated by the force of gravity ($\boldsymbol{F}_g$) and an opposing buoyant force ($\boldsymbol{F}_b$) due to the displaced water. In addition to any translational force moving the centre of mass of the iceberg, where these forces are not colinear they exert a torque, $\boldsymbol{\tau} = \boldsymbol{r} \times \boldsymbol{F}_b$, where $r$ is the vector from $G$ to $B$, which tends to rotate the iceberg.

There are few things one can customize in the code: the shape of the iceberg (defined as a sequence of vertex coordinates), the initial position ($h$) and orientation (call rotate_poly to rotate the iceberg to a given angle), and the friction applied to damp the rotation and translation of the iceberg in water (edit the apply_friction function). I'll tidy up and post the code on GitHub when I get a moment.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

# Plot limits: ideally the figure will not have to rescale during the animation
xmin, xmax = -5, 20
ymin, ymax = -15, 10
width, height = xmax - xmin, ymax - ymin

# The area density of sea water and ice.
rho_water, rho_ice = 1.027, 0.9
# Acceleration due to gravity, m.s-2
g = 9.81

def canonicalize_poly(xy):
"""Shift the (N+1,2) array of coordinates xy to start at minimum y."""
#if get_area(xy) > 0:
#    xy = xy[::-1]
idx_ymin = xy[:,1].argmin()
xy = np.roll(xy[:-1], -idx_ymin, axis=0)
return np.vstack((xy, xy[0]))

def rotate_poly(xy, theta):
"""Rotate the (N+1,2) array of coordinates xy by angle theta about (0,0).

The rotation angle, theta, is in radians.

"""
s, c = np.sin(theta), np.cos(theta)
R = np.array(((c, -s), (s, c)))
xyp = (R @ xy.T).T
return canonicalize_poly(xyp)

def get_area(xy):
"""Return the area of the polygon xy.

xy is a (N+1,0) NumPy array defining the N polygon vertices, but repeating
the first vertex as its last element. The "shoelace algorithm" is used.

"""

x, y = xy.T
return np.sum(x[:-1]*y[1:] - x[1:]*y[:-1]) / 2

def get_cofm(xy, A=None):
"""Return the centre of mass of the polygon xy.

xy is a (N+1,0) NumPy array defining the N polygon vertices, but repeating
the first vertex as its last element. If the polygon area is not passed
in as A it is calculated. The polygon must have uniform density.

"""

if A is None:
A = get_area(xy)
x, y = xy.T
Cx = np.sum((x[:-1] + x[1:]) * (x[:-1]*y[1:] - x[1:]*y[:-1])) / 6 / A
Cy = np.sum((y[:-1] + y[1:]) * (x[:-1]*y[1:] - x[1:]*y[:-1])) / 6 / A
return np.array((Cx, Cy))

def get_moi(xy, rho):
"""Return the moment of inertia of the polygon xy with density rho.

xy is a (N+1,0) NumPy array defining the N polygon vertices, but repeating
the first vertex as its last element.

"""

x, y = xy.T
Ix = rho * np.abs(np.sum((x[:-1]*y[1:] - x[1:]*y[:-1]) *
(y[:-1]**2 + y[:-1]*y[1:] + y[1:]**2)) / 12)
Iy = rho * np.abs(np.sum((x[:-1]*y[1:] - x[1:]*y[:-1]) *
(x[:-1]**2 + x[:-1]*x[1:] + x[1:]**2)) / 12)
# Perpendicular axis theorem.
Iz = Ix + Iy
return Ix, Iy, Iz

def get_zero_crossing(pts):
"""Return the coordinates of the zero-crossing in pts.

pts is a pair of (x, y) points, assumed to be joined by a straight line
segment. This function returns the coordinates (x,0) at which this line
crosses the y-axis (corresponding to sea-level in our model).

"""
P0, P1 = pts
x0, y0 = P0
x1, y1 = P1
if (x1-x0) == 0:
return x1, 0
m = (y1-y0)/(x1-x0)
c = y1 - m*x1
return -c/m, 0

def get_displaced_water_poly(iceberg, submerged=None):
"""Get the polygon for the submerged portion of the iceberg.

iceberg is a (N+1,2) array of coordinates corresponding to the iceberg's
vertexes in its current position and orientation (the first vertex is
repeated at the end of the array);
submerged is a boolean array corresponding to the vertexes which are
under water (<0); if not provided it is calculated.

"""

if submerged is None:
submerged = (iceberg[:,1] <= 0)
nsubmerged = sum(submerged)

# Partially-submerged iceberg: find where it enters the sea, i.e. which
# edges cross zero. zc_idx holds the indexes of the vertices *before*
# each zero-crossing edge.
diff = np.diff(submerged)
zc_idx = np.where(diff)[0]
# Interpolate to find the coordinates of the zero crossing.
ncrossings = len(zc_idx)
# We're going to build a polygon for the shape of the displaced water,
# i.e. the submerged part of the iceberg.
displaced_water = np.empty((nsubmerged + ncrossings, 2))
# Loop over the points *before* each crossing in pairs. NB if the
# iceberg is partially submerged, len(zc_idx) is guaranteed to be even.
assert not ncrossings % 2
i = j  = 0
for idx1, idx2 in zip(zc_idx[0::2], zc_idx[1::2]):
# All the submerged vertices up to the upwards crossing.
displaced_water[j:j+idx1-i+1] = iceberg[i:idx1+1]
# Work out where the crossing vertex should be and add it.
c = get_zero_crossing(iceberg[idx1:idx1+2])
j += idx1 - i + 1
displaced_water[j] = c
j += 1

# Now the downward crossing: all the unsubmerged vertices are
# skipped, and an extra vertex at sea level is added.
c = get_zero_crossing(iceberg[idx2:idx2+2])
displaced_water[j] = c
j += 1
i = idx2 + 1
# Copy across any remaining submerged vertexes to displaced_water.
displaced_water[j:] = iceberg[i:]
return displaced_water

def apply_friction(omega, dh):
"""Apply frictional forces to the angular and linear velocities."""

# Hard friction: angular and linear velocities are immediately quenched.
# after movement.
# return 0, np.array((0,0))

# Intermediate friction: reduce the velocities by some fraction.
return omega * 0.9, dh * 0.6

# Our two-dimensional iceberg, defined as a polygon.
poly = [
(3,0), (3,3), (1,5), (4,7), (0,12), (1,15), (4,17), (6,14), (7,14), (8,12),
(7,10), (7,7), (5,1)
]

# Repeat the initial vertex at the end, for convenience.
iceberg0 = np.array(poly + [poly[0]])
# Centre the iceberg's local coordinates on its centre of mass.
iceberg0 = iceberg0 - get_cofm(iceberg0)
# We might want to start the iceberg off in some different orientation:
# if so, rotate it here.
iceberg0 = rotate_poly(iceberg0, -0.2)

# Get the (signed) area, mass, and weight of the iceberg.
A = get_area(iceberg0)
M = rho_ice * abs(A)
Fg = np.array((0, -M * g))

# We also need the Iz component of the iceberg's moment of inertia.
_, _, Iz = get_moi(iceberg0, rho_ice)

fig, ax = plt.subplots()
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
# We would prefer equal distances in the x- and y-directions to look the same.
ax.axis('equal')

# The centre of mass starts at this height above sea level.
h = 15
# theta is the turning angle of the iceberg from its initial orientation;
# G is the position of its centre of mass (in world coordinates).
theta, G = 0, np.array((6, h))
# omega = dtheta / dt is the angular velocity; dh is the linear velocity.
omega, dh = 0, np.array((0,0))
# The time step (s): small, but not too small.
dt = 0.1

def update(it):
"""Update the animation for iteration number it."""

global omega, dh, G, theta

print('iteration #{}'.format(it))

# Update iceberg orientation and position.
theta += omega * dt
G = G + dh * dt

# Rotate and translate a copy of the original iceberg into its current.
# position.
iceberg = rotate_poly(iceberg0, theta)
iceberg = iceberg + G

# Which vertices are submerged (have their y-coordinate negative)?
submerged = (iceberg[:,1] <= 0)
nsubmerged = sum(submerged)

iceberg_in_water = True
if nsubmerged in (0, 1):
# The iceberg is in the air above the surface of the sea.
B = None
alpha = 0
iceberg_in_water = False

if iceberg_in_water:
# Apply some frictional forces which damp the motion in water.
omega, dh = apply_friction(omega, dh)
if nsubmerged == len(submerged):
# The iceberg is fully submerged.
displaced_water = iceberg
B = G
else:
displaced_water = get_displaced_water_poly(iceberg, submerged)

# Area of the displaced water and position of the centre of buoyancy.
B = get_cofm(displaced_water)

# Buoyant force due to the displaced water.
Fb = np.array((0, rho_water * abs(Adisplaced) * g))

if B is not None:
# Vector from G to B
r = B - G
tau = np.cross(r, Fb)
alpha = tau / Iz

# Resultant force on the iceberg.
F = Fg + Fb
# Net linear acceleration.
a = F / M

# Now plot the scene for this frame of the animation.
ax.clear()

# The sea! The sea!
sea_patch = plt.Rectangle((xmin, ymin), width, -ymin, fc='#8888ff')

# The iceberg itself, in its current orientation and position.
poly_patch = plt.Polygon(iceberg, fc='#ddddff', ec='k')

if B is not None:
# Draw the submerged part of the iceberg in a different colour.
poly_patch = plt.Polygon(displaced_water, fc='#ffdddd', ec='k')

# Indicate the position of the centre of buoyancy.
bofm_patch = plt.Circle(B, 0.2, fc='b')

# Indicate the position of the centre of mass.
cofm_patch = plt.Circle(G, 0.2, fc='r')

# Update the angular and linear velocities
omega += alpha * dt
dh = dh + a * dt

ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.axis('equal')

ani = animation.FuncAnimation(fig, update, 150, blit=False, interval=100, repeat=True)
plt.show()

Current rating: 5

#### John Sharma 3 years, 6 months ago

This is an excellent and informative piece of analysis. Thank you for sharing this, Christian.

Kind Regards,

John

Current rating: 5

#### christian 3 years, 6 months ago

I'm glad you found it interesting – thank you for your comment.

Current rating: 5

#### John Sharma 3 years, 6 months ago

Welcome, Christian. Observing real-life phenomena, like icebergs, being modelled through efficiently written code has been quite fascinating.

John

Current rating: 5

#### Jonas 2 years, 1 month ago

Nice approach,
you may also be interested in the web-applet by Joshua Tauberer https://joshdata.me/iceberger.html