How to plot the motion of a projectile under the effect of gravity, buoyancy and air resistance?

Question

I am trying to make a plot of a projectile motion of a mass which is under the effect of gravitational, buoyancy and drag force. Basically, I want to show effects of the buoyancy and drag force on flight distance, flight time and velocity change on plotting.

import matplotlib.pyplot as plt
import numpy as np

V_initial = 30 # m/s
theta = np.pi/6 # 30
g = 3.711
m =1
C = 0.47
r = 0.5
S = np.pi*pow(r, 2)
ro_mars = 0.0175
t_flight = 2*(V_initial*np.sin(theta)/g)
t = np.linspace(0, t_flight, 200)

# Drag force
Ft = 0.5*C*S*ro_mars*pow(V_initial, 2)

# Buoyant Force
Fb = ro_mars*g*(4/3*np.pi*pow(r, 3))

x_loc = []
y_loc = []

for time in t:
    x = V_initial*time*np.cos(theta)
        y = V_initial*time*np.sin(theta) - (1/2)*g*pow(time, 2)
    x_loc.append(x)
    y_loc.append(y)

x_vel = []
y_vel = []
for time in t:
    vx = V_initial*np.cos(theta)
    vy = V_initial*np.sin(theta) - g*time
    x_vel.append(vx)
    y_vel.append(vy)


v_ch = [pow(i**2+ii**2, 0.5) for i in x_vel for ii in y_vel]

tau = []
for velocity in v_ch:
    Ft = 0.5*C*S*ro_mars*pow(velocity, 2)
        tau.append(Ft)

buoy = []
for velocity in v_ch:
    Fb = ro_mars*g*(4/3*np.pi*pow(r, 3))
    buoy.append(Fb)

after this point, I couldn't figure out how to plot to a projectile motion under this forces. In other words, I'm trying to compare the projectile motion of the mass under three circumstances

1. Mass only under the effect of gravity
2. Mass under the effect of gravity and air resistance
3. Mass under the effect of gravity, air resistance, and buoyancy

Answer 1

You must calculate each location based on the sum of forces at the given time. For this it is better to start from calculating the net force at any time and using this to calculate the acceleration, velocity and then position. For the following calculations, it is assumed that buoyancy and gravity are constant (which is not true in reality but the effect of their variability is negligible in this case), it is also assumed that the initial position is (0,0) though this can be trivially changed to any initial position.

F_x = tau_x 
F_y = tau_y + bouyancy + gravity

Where tau_x and tau_y are the drag forces in the x and y directions, respectively. The velocities, v_x and v_y, are then given by

v_x = v_x + (F_x / (2 * m)) * dt
v_y = v_y + (F_y / (2 * m)) * dt

So the x and y positions, r_x and r_y, at any time t are given by the summation of

r_x = r_x + v_x * dt
r_y = r_y + v_y * dt

In both cases this must be evaluated from 0 to t for some dt where dt * n = t if n is the number of steps in summation.

r_x = r_x + V_i * np.cos(theta) * dt + (F_x / (2 * m)) * dt**2
r_y = r_y + V_i * np.sin(theta) * dt + (F_y / (2 * m)) * dt**2

The entire calculation can actually be done in two lines,

r_x = r_x + V_i * np.cos(theta) * dt + (tau_x / (2 * m)) * dt**2
r_y = r_y + V_i * np.sin(theta) * dt + ((tau_y + bouyancy + gravity) / (2 * m)) * dt**2

Except that v_x and v_y require updating at every time step. To loop over this and calculate the x and y positions across a range of times you can simply follow the below (edited) example.

The following code includes corrections to prevent negative y positions, since the given value of g is for the surface or Mars I assume this is appropriate - when you hit zero y and try to keep going you may end up with a rapid unscheduled disassembly, as we physicists call it.

Edit

In response to the edited question, the following example has been modified to plot all three cases requested - gravity, gravity plus drag, and gravity plus drag and buoyancy. Plot setup code has also been added

Complete example (edited)

import numpy as np
import matplotlib.pyplot as plt

def projectile(V_initial, theta, bouyancy=True, drag=True):
    g = 9.81
    m = 1
    C = 0.47
    r = 0.5
    S = np.pi*pow(r, 2)
    ro_mars = 0.0175

    time = np.linspace(0, 100, 10000)
    tof = 0.0
    dt = time[1] - time[0]
    bouy = ro_mars*g*(4/3*np.pi*pow(r, 3))
    gravity = -g * m
    V_ix = V_initial * np.cos(theta)
    V_iy = V_initial * np.sin(theta)
    v_x = V_ix
    v_y = V_iy
    r_x = 0.0
    r_y = 0.0
    r_xs = list()
    r_ys = list()
    r_xs.append(r_x)
    r_ys.append(r_y)
    # This gets a bit 'hand-wavy' but as dt -> 0 it approaches the analytical solution.
    # Just make sure you use sufficiently small dt (dt is change in time between steps)
    for t in time:
        F_x = 0.0
        F_y = 0.0
        if (bouyancy == True):
            F_y = F_y + bouy
        if (drag == True):
            F_y = F_y - 0.5*C*S*ro_mars*pow(v_y, 2)
            F_x = F_x - 0.5*C*S*ro_mars*pow(v_x, 2) * np.sign(v_y)
        F_y = F_y + gravity

        r_x = r_x + v_x * dt + (F_x / (2 * m)) * dt**2
        r_y = r_y + v_y * dt + (F_y / (2 * m)) * dt**2
        v_x = v_x + (F_x / m) * dt
        v_y = v_y + (F_y / m) * dt
        if (r_y >= 0.0):
            r_xs.append(r_x)
            r_ys.append(r_y)
        else:
            tof = t
            r_xs.append(r_x)
            r_ys.append(r_y)
            break

    return r_xs, r_ys, tof

v = 30
theta = np.pi/4

fig = plt.figure(figsize=(8,4), dpi=300)
r_xs, r_ys, tof = projectile(v, theta, True, True)
plt.plot(r_xs, r_ys, 'g:', label="Gravity, Buoyancy, and Drag")
r_xs, r_ys, tof = projectile(v, theta, False, True)
plt.plot(r_xs, r_ys, 'b:', label="Gravity and Drag")
r_xs, r_ys, tof = projectile(v, theta, False, False)
plt.plot(r_xs, r_ys, 'k:', label="Gravity")
plt.title("Trajectory", fontsize=14)
plt.xlabel("Displacement in x-direction (m)")
plt.ylabel("Displacement in y-direction (m)")
plt.ylim(bottom=0.0)
plt.legend()
plt.show()

Note that this preserves and returns the time-of-flight in the variable tof.

Answer 2

Using vector notation, and odeint.

import numpy as np
from scipy.integrate import odeint
import scipy.constants as SPC
import matplotlib.pyplot as plt

V_initial = 30 # m/s
theta = np.pi/6 # 30
g = 3.711
m = 1 # I assume this is your mass
C = 0.47
r = 0.5
ro_mars = 0.0175

t_flight = 2*(V_initial*np.sin(theta)/g)
t = np.linspace(0, t_flight, 200)

pos0 = [0, 0]
v0 = [np.cos(theta) * V_initial, np.sin(theta)  * V_initial]

def f(vector, t, C, r, ro_mars, apply_bouyancy=True, apply_resistance=True):
    x, y, x_prime, y_prime = vector

    # volume and surface
    V = np.pi * 4/3 * r**3
    S = np.pi*pow(r, 2)

    # net weight bouyancy
    if apply_bouyancy:
        Fb = (ro_mars * V - m) * g *np.array([0,1])
    else:
        Fb = -m  * g * np.array([0,1])

    # velocity vector
    v = np.array([x_prime, y_prime])

    # drag force - corrected to be updated based on current velocity
#    Ft = -0.5*C*S*ro_mars*pow(V_initial, 2)
    if apply_resistance:
        Ft = -0.5*C*S*ro_mars* v *np.linalg.norm(v)
    else:
        Ft = np.array([0, 0])

    # resulting acceleration
    x_prime2, y_prime2 = (Fb + Ft) / m

    return x_prime, y_prime, x_prime2, y_prime2

sol = odeint(f, pos0 + v0 , t, args=(C, r, ro_mars))
plt.plot(sol[:,0], sol[:, 1], 'g', label='tray')
plt.legend(loc='best')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

Note that I corrected your drag force to use the actual (not initial) velocity, I do not know if that was your mistake or it was on purpose.

Also please check the documentation for odeint to understand better how to turn a second order ODE (like the one in your problem) to a first order vector ODE.

To remove air resistance or bouyancy, set apply_bouyancy and apply_resistance to True or False by adding them to the args=(...)