∮ Dynamics of a Body in The Ocean

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.

∮ Non-inertial Reference Frame

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.

∮ Earth as a Non-inertial Frame of Reference

Centrifugal Force

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.

Coriolis Force

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}

Azimuthal Force

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}

Translational Force

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}

Effective Gravity

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}

Geocentric Radius

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.

Buoy oh Buoy

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}

Equations of motion

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.

    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)
        

Trajectory of Object in The Ocean

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} $$


Trajectory of Object in The Ocean

We see that the influence of these fictitious forces are negligible when the existence of much stronger external forces are present.