The question is simple: what are the equations of motion for an object sitting at some point in the ocean? I will be attacking this with a high level of detail in some regards, while ignoring some of the more tricky bits.
Let us go through and analyze what forces will be acting on our object, one at a time.
The first issue to arise in our plight to study objects in the ocean is dealing with this hunk of a non-inertial reference frame we call Earth. The fact that our reference frame is accelerating and rotating is problematic because it introduces the so-called 'fictitious forces'. That is, it appears that forces are acting on objects, but it is just an illusion due to the acceleration of your coordinate system. A prime example is the Coriolis 'force': If a sniper shoots a bullet due North from the equator, it appears that the bullet accelerates to the West. However, this is not due to any force acting on the bullet, it is due to the earth rotating underneath the bullet. Let us explore how these forces arise mathematically.
Consider an arbitrary vector in a moving reference frame
$$
\vec{A} = A_x\hat{x} + A_y \hat{y} + A_z\hat{z}
$$
where our unit vectors change direction over time. Then via the product rule, the time derivative of our vector is
$$
\frac{d\vec{A}}{dt}=\left(\frac{dA_x}{dt}\hat{x} + \frac{dA_y}{dt}\hat{y} + \frac{dA_z}{dt}\hat{z} \right) + \left(A_x\frac{d\hat{x}}{dt} + A_y\frac{d\hat{y}}{dt} + A_z\frac{d\hat{z}}{dt}\right)
$$
The first term on the right hand side is the rate of change of $\vec{A}$ as seen from the axes of our accelerated frame. This is a special quantity, and we denote it as
$$
\frac{\delta \vec{A}}{\delta t} = \frac{dA_x}{dt}\hat{x} + \frac{dA_y}{dt}\hat{y} + \frac{dA_z}{dt}\hat{z}
$$
The final term on the right hand side of our time derivative of $\vec{A}$ is due to the rotation of our axes. Since our axes are rigid, we can say
$$
\frac{d\hat{x}}{dt} = \vec{\omega}\times \hat{x} \qquad \text{and} \qquad \frac{d\hat{y}}{dt} = \vec{\omega}\times \hat{y} \qquad \text{and} \qquad \frac{d\hat{z}}{dt} = \vec{\omega}\times \hat{z}
$$
Where $\vec{\omega}$ is the angular velocity of our rotating coordinate system. So thanks to the cross product operation possesing the distributive property, we can write the time derivative of our arbitrary vector as
$$
\frac{d\vec{A}}{dt} = \frac{\delta \vec{A}}{\delta t} + \vec{\omega}\times \left(A_x\hat{x} + A_y \hat{y} + A_z \hat{z} \right)
$$
\begin{equation}
\frac{d\vec{A}}{dt} = \frac{\delta \vec{A}}{\delta t} + \vec{\omega}\times\vec{A}
\end{equation}
Now say that instead of working with an arbitrary vector, we consider the position of a particle,
$$
\frac{d\vec{r}}{dt} = \frac{\delta \vec{r}}{\delta t} + \vec{\omega}\times \vec{r}
$$
To determine the acceleration of the particle, let's take the time derivative again,
$$
\frac{d^2 \vec{r}}{dt^2}=\frac{d}{dt}\left(\frac{\delta \vec{r}}{dt} + \vec{\omega}\times\vec{r}\right)
$$
Using our results from equation (1), we can rewrite this as
$$
\frac{d^2 \vec{r}}{dt^2}= \frac{\delta}{\delta t}\left( \frac{\delta \vec{r}}{\delta t} + \vec{\omega}\times \vec{r} \right) + \vec{\omega}\times\left( \frac{\delta\vec{r}}{\delta t} + \vec{\omega}\times\vec{r} \right)
$$
Both our derivative operator and the cross product are distributive, so we can rewrite this as
\begin{equation}
\frac{d^2 \vec{r}}{dt^2}= \frac{\delta ^2 \vec{r}}{\delta t^2} + \frac{\delta \vec{\omega}}{\delta t}\times \vec{r} + 2\vec{\omega}\times\frac{\delta \vec{r}}{\delta t} + \vec{\omega}\times(\vec{\omega}\times\vec{r})
\end{equation}
Now we can relate the position of the particle in our accelerated frame to the inertial frame by introducing the vector $\vec{R}$ that stems from the origin of interial frame to the origin of the accelerated frame. We let $\vec{r}_1$ be the position of our particle from our intertial frame of reference:
$$
\vec{r}_1 = \vec{r} + \vec{R}
$$
$$
\frac{d^2\vec{r}_1}{dt^2} = \frac{d^2\vec{r}}{dt^2}+ \frac{d^2\vec{R}}{dt^2}
$$
Using equation (2), this might be rewritten as
$$
\frac{d^2\vec{r}_1}{dt^2} =\frac{\delta ^2 \vec{r}}{\delta t^2} +\dot{\vec{\omega}}\times \vec{r} + 2\vec{\omega}\times\vec{v} + \vec{\omega}\times(\vec{\omega}\times\vec{r})
$$
Where $\dot{\vec{\omega}} = \delta \vec{\omega}/\delta t$ is the angular acceleration of our non-inertial reference frame, and $\vec{v}=\delta \vec{r}/\delta t$ is the velocity of the particle as seen in the non-inertial frame.
Finally, we are able to get our equation of motion for the particle in ou non-inertial frame:
\begin{equation}
m\frac{\delta^2 \vec{r}}{\delta t^2} = \vec{F} - m\left[\vec{\omega}\times(\vec{\omega}\times\vec{r}) + 2\vec{\omega}\times\vec{v} + \dot{\vec{\omega}}\times\vec{r} + \frac{d^2\vec{R}}{dt^2} \right]
\end{equation}
Where $\vec{F}$ is the net force from the inertial frame of reference, and every term in the square brackets is a "fictitious force" that arises in the accelerated frame.
The centrifugal force arises due to the rotation of the earth.
$$
\vec{F}_{cent}= -m\vec{\omega}\times(\vec{\omega}\times\vec{r})
$$
It is true that
$$
\vec{\omega}\cdot\vec{F}_{cent}=0
$$
Hence, this alway acts perpendicularly to the axis of rotation. This is a measurable effect that causes people to weigh less at the equator than they do at the North and South poles.
The angular speed of earth can be found quite easily, since we know that it spins $2 \pi$radians per day:
$$
\omega = 2 \pi\left(\frac{\text{radians}}{\text{day}}\right)\left(\frac{1\text{day}}{24\text{hours}}\right)\left(\frac{1\text{hour}}{3600\text{s}}\right)=7.27\times 10^{-3}\frac{\text{rad}}{\text{s}}
$$
Our coordinate system of choice will have its origin at the surface of the earth, with the $(+)\hat{x}$ direction being East, the $(+)\hat{y}$ direction being North, and the $(+)\hat{z}$ direction being normal to the surface of the earth. This means that the angular velocity of earth will be
\begin{equation}
\vec{\omega} = 0\hat{x} + \omega\sin\theta\hat{y} + \omega\cos\theta\hat{z}
\end{equation}
where $\theta$ is the colatitude angle, measure from the axis of rotation down our y-direction.
The Coriolis force is the apparent deflection of an object due to its velocity in a rotating frame. \begin{equation} \vec{F}_{cori} = -2m\vec{\omega}\times\vec{v} \end{equation}
The azimuthal force arises from the angular acceleration of a coordinate system. Technically, the earth's angular speed is slowing down; this is however, completely negligible for our case. \begin{equation} \vec{F}_{azi} = -m\dot{\vec{\omega}}\times\vec{r} \end{equation}
The translational force arises from the translational acceleration of a non-inertial coordinate system with respect to an inertial coordinate system. Let both of these reference frames have their origins are the center of the earth. Then this force is zero.
But our reference frame is at the surface of the earth, and the origin is moving relative to our non-inertial frame at the center of the earth! This contribution should be negligible.
\begin{equation}
\vec{F}_{trans} = \frac{d^2\vec{R}}{dt^2} = 0
\end{equation}
So, equation (3) now simplifies to
\begin{equation}
m\frac{\delta^2 \vec{r}}{\delta t^2} = \vec{F} - m\left[\vec{\omega}\times(\vec{\omega}\times\vec{r}) + 2\vec{\omega}\times\vec{v} \right]
\end{equation}
If the earth were a perfect sphere of uniform density, then acceleration due to gravity would be the same everywhere on the planet (at equal altitudes). However, due to the centrifugal force previously discussed, the earth is actually a spheroid -- bulging out near the equator. This makes acceleration due to gravity a function of your latitude. We will use the international gravity formula
\begin{equation}
g = g_{eq}\left( 1 + A\sin^2\theta_{lat} -B\sin^2(2\theta_{lat}) \right)
\end{equation}
where $A=0.0053024$ and $B=0.0000058$ are constants, and $g_{eq}=9.780327\frac{\text{m}}{\text{s}^2}$ is acceleration due to gravity at the equation. (Numbers gotten from wikipedia ;).
We combine the contributions of the gravitational force and the centrifugal force to get a quantity called effective gravity.
\begin{equation}
\vec{F}_{grav,eff}=\vec{F}_{grav} - \vec{F}_{cent}
\end{equation}
Since the earth is a spheroid, that means the 'radius', or your distance from the center of the earth, is dependent on your latitude. We use the following equation found on wikipedia. \begin{equation} R_E = \left[\frac{\left(r_{eq}^4\cos^2\theta_{lat} + r_p^4\sin^2\theta_{lat}\right)} { \left(r_{eq}^2\cos^2\theta_{lat} + r_p^2\sin^2\theta_{lat}\right) } \right] \end{equation} where $r_{eq}=6.378137\times 10^6\text{m}$ and $r_p=6.356752\times 10^6 \text{m}$ are the radius at the earth and at the poles, respectively.
We of course need to account for the buoyant force that the ocean water exerts onto our object. It is known that the buoyant force is equal to the weight of water that is displaced by the submerged object. This means,
$$
\vec{F}_{buoy} = \rho V_{sub} g_{eff}
$$
where $\rho$ is the density of water, $V_{sub}$ is the volume of the object that is actually submerged, and $g_{eff}$ is the effective acceleration due to gravity.
We know the density of water, and we have already handled effective acceleration due to gravity, so we need to find the volume of our submerged object, which in this case say is a sphere. Another way to phrase the inquiry of interest is, if we pour liquid into a hollow sphere of radius $R$, what is the volume of liquid within the sphere as a function of its height?
Essentialy we are deriving the volume of a sphere as a function of height. The derivation might be difficult to follow without a visual aid, but such is life.
We will integrate concentric disks from the base of the sphere to its center. Each disk will have a volume element of
$$
dV = \pi(R^2 - (R-h)^2)dh
$$
Integrating over the bottom half of the sphere gives us
$$
V = \pi\int_0^h (R^2 - (R-h)^2)dh = \pi Rh^2 - \frac{\pi}{3}h^3
$$
Of course you now need to integrate the top have of the sphere and add it to the volume of the bottom half. However, things cancel out in such a nice way that we are left with the exact same equation! So,
\begin{equation}
V_{sub} = \pi Rh^2 - \frac{\pi}{3}h^3
\end{equation}
Finally, we have all the pieces needed to find an object's trajectory when it is sitting at some point in the ocean. The equations of motion are as follows
\begin{equation}
\ddot{\vec{r}} = \vec{g}_{eff} - 2\vec{\omega}\times\dot{\vec{r}} + \frac{\vec{F}_{buoy}}{m} + \frac{\vec{F}_{ext}}{m}
\end{equation}
\begin{equation}
\dot{\vec{r}}=\vec{v}
\end{equation}
where $\vec{F}_{ext}$ is the sum of all external forces acting on our object. These could be wind forces, or currents, etc.
Now we will solve these equations of motion using the RK4 method.
I think coding is gross!
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# create a function that does cross products
def cross(a, b):
c_x = a[1]*b[2] - a[2]*b[1]
c_y = a[2]*b[0] - a[0]*b[2]
c_z = a[0]*b[1] - a[1]*b[0]
product = np.array([c_x, c_y, c_z], float)
return(product)
# geocentric radius of earth (function of latitude)
def geocentric_radius(theta):
# theta measured pole->equator ; lat measured equator->pole (lat=pi/2 at pole)
lat = np.pi - theta
top = (r_eq**4)*np.cos(lat)**2 + (r_p**4)*np.sin(lat)**2
bottom = (r_eq**2)*np.cos(lat)**2 + (r_p**2)*np.sin(lat)**2
radius = np.sqrt(top/bottom)
return(radius)
# caluculate rotational angular velocity of earth (w_x, w_y, w_z)
def omega(theta):
w_y = w0*np.sin(theta)
w_z = w0*np.cos(theta)
w_new = np.array([0, w_y, w_z], float)
return(w_new)
# function for gravity as a function of angle (g_x, g_y, g_z)
def gravity(theta):
# theta measured pole->equator ; lat measured equator->pole (lat=pi/2 at pole)
lat = np.pi/2 - theta
grav = g_eq*(1 + A*np.sin(lat)**2 - B*np.sin(2*lat)**2)
grav_ret = np.array([0.0, 0.0, grav], float)
return(grav_ret)
# function for effective gravity (g + centrifugal)
def gravity_effective(g, cent):
grav = g - cent
return(grav)
# buoyancy force
def buoyancy_boi(pos):
depth = max(r_bob - pos[index][2], 0.0) # radius of object minus z-comp, 0 if above sealevel
V_sub = min(np.pi*r_bob*depth**2 - np.pi*depth**3/3, V_bob) # volume submerged
buoyancy = rho*g_effective*V_sub/m
return(buoyancy)
# external force
def force_external(my_time):
return(np.array([0.0, 0.0, 0.0]))
# define a function to calculate the total acceleration; eta = (r, v)
def f(eta):
f_r = eta[1] # slope of position
f_v = (f_ext - corioli - g_effective + buoy) # slope of velocity
return(np.array([f_r, f_v]))
# initial colatitude angle
theta = np.pi/4
# position of object (x, y, z)
r = np.array([[0.0, 0.0, 0.0]], float)
# velocity of object (v_x, v_y, v_z)
v = np.array([[0.0, 0.0, 0.0]], float)
# rotational angular velocity of earth (w_x, w_y, w_z)
w0 = 2*np.pi/(24*3600) # radians/s
w = omega(theta)
# geocentric radius of the earth (function of latitude)
r_eq = 6.378137e6 # meters ; radius at equator
r_p = 6.356752e6 # meters ; radius at poles
# gravity as a function of latitude
g_eq = 9.780327 # m/s^2 ; g at equator
A = 0.0053024
B = 0.0000058
# density of water
rho = 1000 # kg/m^3
# volume of object
r_bob = 0.5 # m ; radius of spherical object
V_bob = (4*np.pi*r_bob**3)/3 # m^3
# other constant and domain
m = 120 # kg
start = 0.0
stop = 10.0
dt = 0.001
steps = int((stop - start)/dt)
domain = np.linspace(start, stop, steps)
# run RK4 algorithm to solve our equations
eta = np.array([r[0], v[0]])
index = 0
time = 0
f_e = force_external(time)
for t in domain:
# radius of earth at our current latitude
rad = geocentric_radius(theta)
# fictitious forces based off our velocity and position
corioli = cross(w, v[index])
centrifug = cross(w, cross(w, r[index]))
# first find acceleration due to gravity as a function of our latitude
# then find the effective gravity by including the centrifugal term
g_effective = gravity_effective(gravity(theta), centrifug)
# buoyancy acceleration
buoy = buoyancy_boi(r)
# exertnal forces (established in our function)
f_ext = force_external(time)
# RK4
k1 = dt*f(eta)
k2 = dt*f(eta + 0.5*k1)
k3 = dt*f(eta + 0.5*k2)
k4 = dt*f(eta + k3)
eta += (k1 + 2*k2 + 2*k3 + k4)/6
# store our updated positions and velocities in our matrices
if len(r) < len(domain):
r = np.vstack([r, eta[0]])
v = np.vstack([v, eta[1]])
f_e = np.vstack([f_e, f_ext])
index += 1
time += dt
# update the colatitude angle
theta += dt*v[index][1]/(rad*np.sin(theta))
# update components of rotational angular speed of earth
w = omega(theta)
else:
break
### now animate ###
# stylize
with plt.rc_context({'axes.edgecolor':'white', 'xtick.color':'white', \
'ytick.color':'white', 'figure.facecolor':'darkslategrey',\
'axes.facecolor':'darkslategrey','axes.labelcolor':'white',\
'axes.titlecolor':'white'}):
fig = plt.figure(figsize=(10,8))
fig.tight_layout()
ax = plt.axes(projection ='3d')
ax.view_init(10, -90) # view angles of 3dplot
#ax.view_init(90, -90) # top view
# plot water level as a plane
x = np.linspace(min(r[:,0]), max(r[:,0]),10)
y = np.linspace(min(r[:,1]), max(r[:,1]),10)
X,Y = np.meshgrid(x,y)
Z= 0 + 0*X + 0*Y
ax.plot_surface(X, Y, Z, alpha=0.3)
# axes labels and rotate tick markings
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.xticks(rotation=30)
ims = [] # images list for animation
index = 0 # what timestep are we on
update_text = False # updating text with strength of f_external
text_update = 0 # to help limit how quickly text updates
while index < len(domain):
# plot the trajectory up to our timestep
curve, = ax.plot3D(r[:index,0], r[:index,1], r[:index,2], 'm')
# plot the center of mass of the body
bob = ax.scatter(r[index,0], r[index,1], r[index,2], color='orange', s=1000)
# occasionally update some text that tells us f_ext
if (text_update%10 == 0) and (update_text == True):
txt = ax.text(0.001, 2, -0.01, s=str(round(f_e[index,1],3))+'N north')
# add our pictures to the ims list
ims.append([bob, curve, txt])
# only show every 20 time steps
index += 20
text_update += 1
ani = animation.ArtistAnimation(fig, ims, interval=100)
writervideo = animation.FFMpegWriter(fps=60)
ani.save('./bobInWater.mp4', writer=writervideo)
The orange blob represents the center of mass of the body. We see that the object bounces up and down as the buoyant force and the force of gravity battle one another. We see a negligible amount of movement in the y-direction (on the order of 40nm North or so), but a more significant amount of movement in the x-direction (a shocking 0.1mm West).
Now let us add a time dependent force acting in the North direction, of the form
$$
\vec{F}_{ext} = 0\hat{x} + A\sin^2 (t) \hat{y} + 0 \hat{z}
$$
We see that the influence of these fictitious forces are negligible when the existence of much stronger external forces are present.