Path: Notes/Day 5: Drag / Drag.ipynb

Views: ²⁰²

Kernel: VPython

ParseError: KaTeX parse error: Undefined control sequence: \require at position 1: \̲r̲e̲q̲u̲i̲r̲e̲{cancel}ParseError: KaTeX parse error: Undefined control sequence: \require at position 1: \̲r̲e̲q̲u̲i̲r̲e̲{enclose}

Drag

PY345

Puzzle of the Day

Drag

air resistance is also known as drag
drag points opposite the direction of motion
drag is velocity-dependent, $\vec{F}^d=-f(v)\hat{v}$
in general $f(v)$ can be complicated, but for low speeds it is approximately true that $\vec{F}^d=(-bv-cv^2)\hat{v},$ where $b$ and $c$ are constants that depend on the size and shape of the object as well as the properties of the medium through which it is traveling and $\hat{v}=\frac{\vec{v}}{|\vec{v}|}$

Practice

A cart is traveling horizontally with negligible friction. Find horizontal position and speed as a function of time for the case $bv\gg cv^2$ .

Solution

Assumptions

the cart has mass $m$
friction and $cv^2$ are negligible
the cart starts at $x(0)=x_0$ and with speed $v(0)=v_0$

Diagrams

Analysis

According to the free body diagram, Newton's Second Law in the x-direction is

$-F{ac}^d=m\ddot{x}$

$-bv=m\ddot{x}$ .

We can solve this by recognizing that $\ddot{x}=\dot{v}$ .

$-bv=m\dot{v}=m\frac{dv}{dt}$

Now we get v's on one side and t's on the other (it doesn't matter what side the constants are on).

$-\frac{b}{m}dt=\frac{dv}{v}$

Integrate both sides.

$-\int\limits_0^t\frac{b}{m}dt'=\int\limits_{v_0}^{v(t)}\frac{dv}{v}$

$-\frac{b}{m}t=\ln\left (\frac{v(t)}{v_0}\right )$

Exponentiate both sides.

$e^{-\frac{b}{m}t}=\frac{v(t)}{v_0}$

$v(t)=v_0e^{-\frac{b}{m}t}$

Integrate this to find $x(t)$ .

$v(t)=\frac{dx}{dt}=v_0e^{-\frac{b}{m}t}\Rightarrow dx=v_0e^{-\frac{b}{m}t}dt$

$\int\limits_{x_0}^{x(t)}dx=\int\limits_0^tv_0e^{-\frac{b}{m}t'}dt'$

$x(t)-x_0=-\frac{m}{b}v_0\left [e^{-\frac{b}{m}t'}\right ]_0^t=-\frac{mv_0}{b}\left (e^{-\frac{b}{m}t}-1\right )$

$x(t)=-\frac{mv_0}{b}\left (e^{-\frac{b}{m}t}-1\right )+x_0$

Check

Let's check units first. Since $bv$ must have units of force, the SI units of $b$ must be $\frac{[N]}{[m/s]}$ . That means that $\frac{b}{m}t$ has units of $\frac{[N]/[m/s]}{[kg]}[s]=\frac{[kg][m]/[s]^2/[m/s]}{[kg]}[s]=\text{no units}$ , which is good since the argument of an exponential can't have units.

The SI units of $\frac{mv_0}{b}$ are $\frac{[kg][m/s]}{[N]/[m/s]}=\frac{[kg][m/s]^2}{[N]}=\frac{[N][m]}{[N]}=[m]$ , which is what we expect since it is the coefficient in front of a unitless function in the $x(t)$ equation.

Now check limits. ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to 0}… and ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to 0}…, which are expected. ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to \i… and ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to \i…. The former is expected since drag will bring the cart to a stop. The latter seems plausible, as it is at least bigger than $x_0$ .

Interpretation

Linear drag causes the velocity to decay exponentially.

In [1]:

from __future__ import division, print_function
from vpython import *

#initial values
x=0
v=10
m=1
b=1
t=0
dt=0.01

#draw objects
scene=canvas(title="horizontal motion with linear drag")
cart=box(pos=vec(x,0,0), length=1, width=0.5, height=0.5)

#create graphs
xgraph=graph(xtitle="t (s)", ytitle="x (m)")
xdata=gcurve(color=color.blue)

#animate motion
while v>=0.1: #continue looping until cart reaches a very small speed (it asymptotically approaches zero)
    rate(1/dt)
    a=-b*v/m #calculate acceleration from Newton's Second Law
    v=v+a*dt #update velocity
    cart.pos.x=cart.pos.x+v*dt #update position
    t=t+dt #increment time
    xdata.plot(t,cart.pos.x) #plot x vs. t

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

Practice

A cart is traveling horizontally with negligible friction. Find horizontal position and speed as a function of time for the case $cv^2\gg bv$ .

Solution

Assumptions

the cart has mass $m$
friction is negligle, but there is air resistance for which $cv^2\gg bv$
the cart starts at $x(0)=x_0$ with initial velocity $v(0)=v_0$

Diagrams

Analysis

According to the free body diagram, Newton's Second Law in the x-direction is

$-F_{ac}^d=m\ddot{x}$

$-cv^2=m\dot{v}=m\frac{dv}{dt}$ .

Move all of the t's to one side and the v's to the other side and integrate.

$-\frac{c}{m}dt=\frac{dv}{v^2}$

$-\frac{c}{m}\int\limits_0^tdt'=\int\limits_{v_0}^{v(t)}\frac{dv}{v^2}$

$-\frac{c}{m}t=\left [-\frac{1}{v}\right ]_{v_0}^{v(t)}=-\frac{1}{v(t)}+\frac{1}{v_0}$

$\frac{1}{v(t)}=\frac{1}{v_0}+\frac{c}{m}t$

$v(t)=\frac{1}{\frac{1}{v_0}+\frac{c}{m}t}=\frac{v_0}{1+\frac{v_0c}{m}t}$

To find $x(t)$ , we integrate.

$v(t)=\frac{dx}{dt}=\frac{v_0}{1+\frac{v_0c}{m}t}\Rightarrow dx=\frac{v_0}{1+\frac{v_0c}{m}t}dt$

If it weren't for the 1+ in the denominator, we would identify this as an integral that results in a natural log. We can make a change of variables though to make it look like that.

Let $u=1+\frac{v_0c}{m}t$ , then $du=\frac{v_0c}{m}dt\Rightarrow dt=\frac{m}{v_0c}du$ . Then our integral becomes

$\int\limits_{x_0}^{x(t)}dx=\int\limits_{u_0}^{u(t)}\frac{v_0}{u}\frac{m}{v_0c}du=\frac{m}{c}\int\limits_{u_0}^{u(t)}\frac{du}{u}$

$x(t)-x(0)=\frac{m}{c}\left [\ln u\right ]_{u_0}^{u(t)}=\frac{m}{c}\ln\left (\frac{u(t)}{u_0}\right )$

Now we put it back in terms of $t$ by using $u(t)=1+\frac{v_0c}{m}t$ and $u(0)=u_0=1+\frac{v_0c}{m}(0)=1$ .

$x(t)=\frac{m}{c}\ln(1+\frac{v_0c}{m}t)+x_0$

Check

First, let's check units. Since $cv^2$ must have units of force, then $c$ must have SI units of $\frac{[N]}{[m/s]^2}=\frac{[kg][m]/[s]^2}{[m]^2/[s]^2}=\frac{[kg]}{[m]}$ . $\frac{v_0c}{m}t$ then has SI units of $\frac{[m/s][kg/m]}{[kg]}[s]=\text{no units}$ , which is expected since the argument of a natural log function can't have units (and also since it is being added to a unitless number). The coefficient $\frac{m}{c}$ in front of the natural log in $x(t)$ has SI units $\frac{[kg}{[kg]/[m]}=[m]$ , which is what we expect.

Now check limits. ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to 0}… and ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to 0}…, which are both expected. ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to \i… and ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to\in…. The former is expected since drag brings it to a stop, but the latter is surprising because it comes to a stop.

Interpretation

Apparently even though the cart constantly slows down, it only ever comes to a stop after an infinite amount of time.

In [2]:

from __future__ import division, print_function
from vpython import *

#initial values
x=0
v=10
m=1
c=1
t=0
dt=0.01

#draw objects
scene=canvas(title="horizontal motion with linear drag")
cart=box(pos=vec(x,0,0), length=1, width=0.5, height=0.5)

#create graphs
xgraph=graph(xtitle="t (s)", ytitle="x (m)")
xdata=gcurve(color=color.blue)

#animate motion
while v>=0.1: #continue looping until cart reaches a very small speed (it asymptotically approaches zero)
    rate(1/dt)
    a=-c*v**2/m #calculate acceleration from Newton's Second Law
    v=v+a*dt #update velocity
    cart.pos.x=cart.pos.x+v*dt #update position
    t=t+dt #increment time
    xdata.plot(t,cart.pos.x) #plot x vs. t

<IPython.core.display.Javascript object>

Practice

A sphere is dropped from rest from a very large tower. Find height and speed as a function of time for the sphere for the case $bv\gg cv^2$ .

Solution

Assumptions

the $cv^2$ term in the drag equation is negligible
the sphere has a mass $m$

Diagrams

Analysis

According to the free body diagram

$F_{as}^d-F_{es}^g=m\ddot{y}$

$-bv-mg=m\ddot{y}=m\frac{dv}{dt}$ .

You may be wondering why there is a negative in front of the $bv$ term since the force points in the upward direction. In fact, the first time I did this problem that's exactly what I did (I discovered the mistake by checking my answer). But recall that $v$ can be either positive or negative. If it is positive, meaning upward, then the direction of the drag force will be downward (which is what we want). If $v$ is negative, meaning downward, then the direction of the drag force will be upward.

Our goal is to get all of the v's on one side and all of the t's on the other side.

$dt=\frac{mdv}{-bv-mg}$

If the $-mg$ weren't there, this would be a relatively simple integral. In order to solve it, we are going to be a bit creative. First, factor out $-b$ from the bottom.

$dt=\frac{mdv}{-b\left (v+\frac{mg}{b}\right )}$

Now let $u=v+\frac{mg}{b}$ , which means that $du=dv$ . This is called a change of variables.

$dt=-\frac{m}{b}\frac{du}{u}$

Now we integrate both sides.

$\int\limits_0^tdt'=-\frac{m}{b}\int\limits_{u_0}^{u(t)}\frac{du}{u}$

$t=-\frac{m}{b}\ln\left (\frac{u(t)}{u_0}\right )$

Exponentiate both sides.

$e^t=e^{-m/b}\left (\frac{u(t)}{u_0}\right )$

$e^{-\frac{b}{m}t}=\frac{u(t)}{u_0}\Rightarrow u(t)=u_0e^{-\frac{b}{m}t}$

Now let's put it back in terms of $v$ instead of $u$ .

$v(t)+\frac{mg}{b}=\left (v_0+\frac{mg}{b}\right )e^{-\frac{b}{m}t}$

$v(t)=\left (v_0+\frac{mg}{b}\right )e^{-\frac{b}{m}t}-\frac{mg}{b}$

To find $y(t)$ , we integrate with respect to $t$ .

$\int\limits_{y_0}^{y(t)}dy=\int\limits_0^t\left [\left (v_0+\frac{mg}{b}\right )e^{-\frac{b}{m}t}-\frac{mg}{b}\right]dt'$

$y(t)-y_0=\left [\left (v_0+\frac{mg}{b}\right )e^{-\frac{b}{m}t'}-\frac{mg}{b}t' \right]_0^t$

$y(t)-y_0=-\frac{m}{b}\left (v_0+\frac{mg}{b}\right )e^{-\frac{b}{m}t}-\frac{mg}{b}t+\frac{m}{b}\left (v_0+\frac{mg}{b}\right)$

$y(t)=\frac{m}{b}\left (v_0+\frac{mg}{b}\right )\left (1-e^{-\frac{b}{m}t}\right )-\frac{mg}{b}t+y_0$

Check

$bv$ must have units of force, which means that $b$ has units of force/(distance/time). Any term involving $\frac{mg}{b}$ , then must have units of $\frac{\text{force}}{\text{force}/(\text{distance}/\text{time})}=\frac{\text{distance}}{\text{time}}$ . $\frac{b}{m}t$ , meanwhile, has units of $\frac{\text{force}/(\text{distance}/\text{time})\times \text{time}}{\text{mass}}$ . Since $F/m=a$ , $\frac{b}{m}t$ has units of $\frac{\text{distance}}{\text{time}^2}\times \frac{\text{time}^2}{\text{distance}}$ , which has no units. This is good, since the argument of an exponential can never have units.

Now let's check the units in $y(t)$ . We have already established that $\frac{mg}{b}$ has units of speed. That means that $\frac{m}{b}$ has units of speed/acceleration. Therefore $\frac{mg}{b}\times\frac{m}{b}$ has units of $\text{speed}\times\frac{\text{speed}}{\text{acceleration}}=\frac{\text{speed}^2}{\text{acceleration}}=\frac{\text{distance}^2/\text{time}^2}{\text{distance}/\text{time}^2}=\text{distance}$ , which is what we expect.

Let's also check $v(0)$ and $y(0)$ .

$v(0)=\left (v_0+\frac{mg}{b}\right )-\frac{mg}{b}=v_0$

$y(0)=\frac{m}{b}\left (v_0+\frac{mg}{b}\right )\left (1-1\right )-0+y_0=y_0$

As time goes on, $e^{-\frac{b}{m}t}$ gets smaller, which means that $v(t)$ starts at zero and increases in the negative direction. $y(t)$ , meanwhile becomes smaller and smaller as time goes on because the coefficient in front of the parentheses is negative. Both of these behaviors are expected.

Interpretation

If $t$ is very large, then ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to \i…, which means that ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to \i…. This is referred to as terminal velocity.

In [1]:

from __future__ import division, print_function
from vpython import *

#height from which ball is dropped
h=100

#draw objects
scene=canvas(center=vec(0,h/2,0))
building=box(pos=vec(0,h/2,0),length=10,width=10,height=h)
ball=sphere(pos=vec(7,h,0),color=color.yellow,radius=1)

#set initial values
m=1 #mass
g=9.8 #gravitational field strength
b=1 #linear drag coefficient
v=0 #initial velocity
t=0 #initial time
dt=0.01 #time step

#create graphs of height vs. time
ygraph=graph(title="y vs. t (blue=numerical with drag, red=no drag, green=exact with drag)", xtitle="t (s)", ytitle="y (m)", ymin=0, ymax=h)
ydata=gcurve(color=color.blue) #height vs. time with drag
nodrag=gcurve(color=color.red) #height vs. time without drag
yexact=gcurve(color=color.green) #exact solution of height vs. time

#create graphs of velocity vs. time
vgraph=graph(title="v vs. t (blue=numerical with drag, green=exact with drag)", xtitle="t (s)", ytitle="v (m/s)")
vdata=gcurve(color=color.blue) #velocity vs. time
vexact=gcurve(color=color.green) #exact solution of velocity vs. time


while ball.pos.y>=0: #continue while ball is above ground level
    rate(1/dt) #make sure the looping rate doesn't exceed real-time
    a=-b*v/m-g #acceleration from Newton's Second Law
    v=v+a*dt #update velocity
    ball.pos.y=ball.pos.y+v*dt #update height
    t=t+dt #increment time
    #plot all the things
    ydata.plot(t,ball.pos.y) 
    nodrag.plot(t,h-0.5*g*t**2)
    yexact.plot(t,-m*g/b*(m/b*exp(-b/m*t)+t-m/b)+h)
    vdata.plot(t,v)
    vexact.plot(t,m*g/b*(exp(-b/m*t)-1))

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

Practice

A sphere is dropped from rest from a very large tower. Find height and speed as a function of time for the sphere for the case $cv^2\gg bv$ .

Solution

Assumptions

the $bv$ term in the drag equation is negligible
the sphere has mass $m$

Diagrams

Analysis

According to the free body diagram

$F_{as}^d-F_{es}^g=m\ddot{y}$

$cv^2-mg=m\dot{v}=m\frac{dv}{dt}$

Notice that in the linear case there was a negative sign in front of the drag term, but this time there isn't. What gives? The problem is that the sign of $cv^2$ doesn't change when the direction of the velocity changes because the velocity gets squared. So here I am making it positive since I know that the drag force will always be pointing upward. (If this difference in the way you use linear and quadratic drag distresses you, you can multiply the drag term by $\frac{v}{|v|}$ and then always use the negative.)

Let's try using the same tactic that we did in the linear case. Start by factoring out $c$ .

$c\left (v^2-\frac{mg}{c}\right )=m\frac{dv}{dt}$

Then we move all of the $t$ terms to one side and all of the $v$ terms to the other side (it doesn't matter what side the constants are on).

$\frac{c}{m}dt=\frac{dv}{v^2-\frac{mg}{c}}$

Then integrate.

$\int\limits_0^t\frac{c}{m}dt'=\int\limits_0^{v(t)}\frac{dv}{v^2-\frac{mg}{c}}$

The last trick we used in the linear case was to institute a change of variables. It doesn't work out quite the same way here. Let me show you. Let $u=v^2-\frac{mg}{c}$ , then $du=2vdv\Rightarrow dv=\frac{1}{2v}du$ .

$\int\limits_0^t\frac{c}{m}dt'=\int\limits_{u_0}^{u(t)}\frac{\frac{1}{2v}du}{u}$

Notice that there is still a $v$ left over. So if not a change of variables, what's our next option? How about looking up the integral in an integral table. Here is what I find.

$\int\frac{dx}{x^2-a^2}=-\frac{1}{a}\coth^{-1}\left (\frac{x}{a}\right )$ if $x^2>a^2$

In our case, $x=v$ and $a=\sqrt{\frac{mg}{c}}$ , but it is not true that $x^2>a^2$ since $v$ starts at zero. A related integral from the integral table is

$\int\frac{dx}{a^2-x^2}=\frac{1}{a}\tanh^{-1}\left (\frac{x}{a}\right )$ if $x^2<a^2$ . We can make our integral look like this one by multiplying both sides of our equation by $-1$ .

$-\int\limits_0^t\frac{c}{m}dt'=-\int\limits_0^{v(t)}\frac{dv}{v^2-\frac{mg}{c}}=\int\limits_0^{v(t)}\frac{dv}{\frac{mg}{c}-v^2}$

For brevity, let $v_t=\sqrt{\frac{mg}{c}}$ .

$-\frac{c}{m}t=\left [\frac{1}{v_t}\tanh^{-1}\left (\frac{v}{v_t}\right )\right ]_0^{v(t)}$

ParseError: KaTeX parse error: Undefined control sequence: \cancelto at position 60: …}{v_t}\right )-\̲c̲a̲n̲c̲e̲l̲t̲o̲{0}{\tanh^{-1}(…

$v(t)=v_t\tanh\left (-\frac{cv_t}{m}t\right )=-v_t\tanh\left (\frac{cv_t}{m}t\right )$

We integrate to find $y(t)$ .

$\int\limits_{y_0}^{y(t)}dy=-\int\limits_0^tv_t\tanh\left (\frac{cv_t}{m}t'\right )dt'$

Again, I just looked up this integral in a table.

$y(t)=-v_t\left (\frac{1}{\frac{cv_t}{m}}\right )\left [\ln\cosh\left (\frac{cv_t}{m}t'\right )\right ]_0^t+y_0$

$y(t)=-\frac{m}{c}\ln\left [\cosh\left (\frac{cv_t}{m}t\right )-0\right ] +y_0=-\frac{m}{c}\ln\left [\cosh\left (\frac{cv_t}{m}t\right )\right ] +y_0$

Check

$\frac{cv_t}{m}t=\frac{c}{m}\sqrt{\frac{mg}{c}}=\sqrt{\frac{mgc^2}{m^2c}}t=\sqrt{\frac{gc}{m}}t$ . Since the drag force is $cv^2$ , the SI units of $c$ are $c=\frac{[N]}{[m/s]^2}=\frac{[kg\cdot m/s^2]}{[m^2/s^2]}=\frac{[kg]}{[m]}$ . Therefore $\frac{cv_t}{m}t=\sqrt{\frac{gc}{m}}t$ has units of $\sqrt{\frac{[m/s^2][kg/m]}{[kg]}}[s]=\sqrt{\frac{1}{[s]^2}}[s]=\text{no units}$ , which is what we expect since the argument of a trig function can't have units.

The units of $v_t=\sqrt{\frac{mg}{c}}$ are $\sqrt{\frac{[kg][m/s^2]}{[kg/m]}}=\sqrt{\frac{[m]^2}{[s]^2}}=\frac{[m]}{[s]}$ , which is what we expected since that is the coefficient out front of a trig term in the $v(t)$ equation.

In the $y(t)$ equation, the coefficient out front is $\frac{m}{c}$ . This has units of $\frac{[kg]}{[kg/m]}=[m]$ , which is what we expect.

Now let's check limits. ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to 0}…, which is good since the sphere is dropped from rest and that makes $v(0)=0$ . ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to\in…, which is also good since that means that ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to\in…. That is, the velocity approaches a constant, which is consistent with the idea that object's experiencing drag have a terminal velocity.

ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to 0}…, and $\ln(1)=0$ , so ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to 0}…, which is what we expect. Meanwhile, ParseError: KaTeX parse error: Expected group after '_' at position 5: \lim_̲\limits{t\to \i…, which makes sense mathematically since it will just keep falling (we didn't incorporate the ground in any of the analysis.

Interpretation

While the functions are more complicated in the case of quadratic drag vs. linear drag, it still exhibits a terminal velocity.

In [1]:

from __future__ import division, print_function
from vpython import *

h=100 #height from which ball is dropped

#draw objects
scene=canvas(center=vec(0,h/2,0))
building=box(pos=vec(0,h/2,0),length=10,width=10,height=h)
ball=sphere(pos=vec(7,h,0),color=color.yellow,radius=1)

#initial values
m=1 #mass
g=9.8 #gravitational field strength
c=1 #quadratic drag coefficient
v=0 #initial velocity
t=0 #initial time
dt=0.01 #time step

#create height vs. time graphs
ygraph=graph(title="y vs. t (blue=numerical with drag, red=no drag, green=exact with drag)", xtitle="t (s)", ytitle="y (m)", ymin=0, ymax=h)
ydata=gcurve(color=color.blue) #height vs. time with drag
nodrag=gcurve(color=color.red) #height vs. time without drag
yexact=gcurve(color=color.green) #exact solution of height vs. time

#create velocity vs. time graphs
vgraph=graph(title="v vs. t (blue=numerical with drag, green=exact with drag)", xtitle="t (s)", ytitle="v (m/s)")
vdata=gcurve(color=color.blue) #velocity vs. time
vexact=gcurve(color=color.green) #exact solution of velocity vs. time

while ball.pos.y>=0: #continue while ball is above ground
    rate(1/dt) #make sure looping rate does not exceed real-time
    a=c*v**2/m-g #calculate acceleration from Newton's Second Law
    v=v+a*dt #update velocity
    ball.pos.y=ball.pos.y+v*dt #update height
    t=t+dt #increment time
    #plot all the things
    ydata.plot(t,ball.pos.y)
    nodrag.plot(t,h-0.5*g*t**2)
    yexact.plot(t,-m/c*log(cosh(c*sqrt(m*g/c)*t/m))+h)
    vdata.plot(t,v)
    vexact.plot(t,-sqrt(m*g/c)*tanh(c*sqrt(m*g/c)*t/m))

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

Practice

A baseball is thrown upward at an angle. Find position and velocity as a function of time. Assume linear air resistance.

Solution

Assumptions

the ball starts at $x(0)=x_0$ and $y(0)=y_0$ with velocity $v_x(0)=v_{0x}$ and $v_y(0)=v_{0y}$
there is air resistance, where $bv\gg cv^2$
the ball has mass $m$

Diagrams

Analysis

Let's start by writing Newton's Second Law in full vector form.

$\vec{F}_{ab}^d+\vec{F}_{eb}^g=m\ddot{\vec{r}}$

$-bv\hat{v}-mg\hat{y}=m\ddot{x}\hat{x}+m\ddot{y}\hat{y}$

Recall that the way you get a unit vector is to divide a vector by its own magnitude. In this case, $\hat{v}=\frac{\vec{v}}{v}=\frac{v_x\hat{x}+v_y\hat{y}}{v}$ .

$-b\cancel{v}\left (\frac{v_x\hat{x}+v_y\hat{y}}{\cancel{v}}\right )-mg\hat{y}=m\ddot{x}\hat{x}+m\ddot{y}\hat{y}$

We can now separate it into two separate equation. First, the x-direction:

$-bv_x=m\ddot{x}$

But we have already solved this equation. It is the same as the horizontally moving cart with linear drag. Rather than reworking it here, I quote the solution.

$v_x(t)=v_{0x}e^{-\frac{b}{m}t}$ and $x(t)=-\frac{mv_{0x}}{b}\left (e^{-\frac{b}{m}t}-1\right )+x_0$

For the y-direction:

$-bv_y-mg=m\ddot{y}$

We have also already done this problem. It is the same as the equation we derived for the falling sphere.

$v_y(t)=\left (v_{0y}+\frac{mg}{b}\right )e^{-\frac{b}{m}t}-\frac{mg}{b}$ and $y(t)=\frac{m}{b}\left (v_0+\frac{mg}{b}\right )\left (1-e^{-\frac{b}{m}t}\right )-\frac{mg}{b}t+y_0$

Check

We checked these in previous problems.

Interpretation

For linear drag, projectile motion reduces to two one-dimensional problems, which means that the horizontal and vertical motions are independent of one another.

As with the falling objects we discussed before, our solution no longer applies once the object hits the ground.

In [24]:

from __future__ import division, print_function
from vpython import *

#initial conditions
x0=0 #initial x-coordinate of ball
y0=0 #initial y-coordinate of ball
v0x=10 #initial x-component of velocity
v0y=10 #initial y-component of velocity
g=9.8 #gravitational field strength
b=0.1 #drag coefficient
t=0 #initial time
dt=0.01 #time step

#draw objects
scene=canvas(title="projectile motion with linear drag", center=vec(5,2,0))
ball=sphere(pos=vec(x0,y0,0),v=vec(v0x,v0y,0),m=1, color=color.white, radius=0.1, make_trail=True)
ground=box(pos=vec(10,-0.1,0),length=20, width=10,height=0.1,color=color.green)

#create graphs
xygraph=graph(title="trajectory of ball (blue=numerical, green=exact, red=no drag)", xtitle="x (m)", ytitle="y (m)")
xydata=gcurve(color=color.blue)
xynodrag=gcurve(color=color.red)
xyexact=gcurve(color=color.green)

#define gravitational force
Fg=vec(0,-ball.m*g,0)


while ball.pos.y>=0: #continue while ball is above ground level
    rate(1/dt) #make sure looping rate does not exceed real time
    a=(-b*mag(ball.v)*hat(ball.v)+Fg)/ball.m #calculate acceleration from Newton's Second Law
    ball.v=ball.v+a*dt #update velocity
    ball.pos=ball.pos+ball.v*dt #update position
    t=t+dt #increment time
    #plot all the things
    xydata.plot(ball.pos.x, ball.pos.y)
    xyexact.plot(-ball.m*v0x/b*(exp(-b*t/ball.m)-1)+x0,ball.m/b*(v0y+ball.m*g/b)*(1-exp(-b*t/ball.m))-ball.m*g/b*t+y0)
    xynodrag.plot(x0+v0x*t,y0+v0y*t-0.5*g*t**2)

<IPython.core.display.Javascript object>

Practice

A baseball is thrown upward at an angle. Find position and velocity as a function of time. Assume quadratic air resistance.

Solution

Assumptions

the ball starts at $x(0)=x_0$ and $y(0)=y_0$ with velocity $v_x(0)=v_{0x}$ and $v_y(0)=v_{0y}$
there is air resistance, where $cv^2\gg bv$
the ball has mass $m$

Diagrams

Analysis

According to the free body diagram:

$\vec{F}_{ab}^d+\vec{F}_{eb}^g=m\ddot{\vec{r}}$

$-cv^2\hat{v}-mg\hat{y}=m\ddot{x}\hat{x}+m\ddot{y}\hat{y}$

$-c\left (v_x^2+v_y^2\right)\left (\frac{v_x\hat{x}+v_y\hat{y}}{\sqrt{v_x^2+v_y^2}}\right )-mg\hat{y}=m\ddot{x}\hat{x}+m\ddot{y}\hat{y}$

$-c\sqrt{v_x^2+v_y^2}\left (v_x\hat{x}+v_y\hat{y}\right )-mg\hat{y}=m\ddot{x}\hat{x}+m\ddot{y}\hat{y}$

Now we separate into two separate equations:

$-c\sqrt{v_x^2+v_y^2}v_x=m\ddot{x}=m\dot{v_x}$

and

$-c\sqrt{v_x^2+v_y^2}v_y-mg=m\ddot{y}=m\dot{v_y}$ .

Our usual method of putting all of the v's on one side and putting all of the t's on the other side won't work here. The reason is because $v_x$ and $v_y$ both appear in both equations. The problem with that is that if you try to integrate over $dv_x$ , for example, the value of $v_y$ is not constant over that integration. We say that these two differential equations are coupled. Coupled equations can be solved sometimes, but in this particular scenario there is no analytic solution. The solution must be found numerically.

Check

Even though we don't have solutions for $x(t)$ and $y(t)$ , there are still some things we can do to check. For example, we can suppose that $v_y(t)=0$ . In that case, the x-equation becomes

$-cv_x^2=m\ddot{x}$ ,

which is the same equation we had for the cart moving horizontally with quadratic drag.

What if $v_x(t)=0$ ? In that case,

$-c(\pm v_y)v_y-mg=m\ddot{y}$ .

In this case I was careful to use $\pm$ on the square root for the reasons that we mentioned before with the falling sphere with quadratic drag. On the way up, we use the positive solution (which makes the drag force point down), and on the way down we use the negative solution (which makes the drag force point up).

Interpretation

In [27]:

from __future__ import division, print_function
from vpython import *

#initial conditions
x0=0 #initial x-coordinate of ball
y0=0 #initial y-coordinate of ball
v0x=10 #initial x-component of velocity
v0y=10 #initial y-component of velocity
g=9.8 #gravitational field strength
c=0.01 #drag coefficient
t=0 #initial time
dt=0.01 #time step

#draw objects
scene=canvas(title="projectile motion with quadratic drag", center=vec(5,2,0))
ball=sphere(pos=vec(x0,y0,0),v=vec(v0x,v0y,0),m=1, color=color.white, radius=0.1, make_trail=True)
ground=box(pos=vec(10,-0.1,0),length=20, width=10,height=0.1,color=color.green)

#create graphs
xygraph=graph(title="trajectory of ball (blue=numerical, red=no drag)", xtitle="x (m)", ytitle="y (m)")
xydata=gcurve(color=color.blue)
xynodrag=gcurve(color=color.red)

#define gravitational force
Fg=vec(0,-ball.m*g,0)


while ball.pos.y>=0: #continue while ball is above ground level
    rate(1/dt) #make sure looping rate does not exceed real time
    a=(-c*mag(ball.v)**2*hat(ball.v)+Fg)/ball.m #calculate acceleration from Newton's Second Law
    ball.v=ball.v+a*dt #update velocity
    ball.pos=ball.pos+ball.v*dt #update position
    t=t+dt #increment time
    #plot all the things
    xydata.plot(ball.pos.x, ball.pos.y)
    xynodrag.plot(x0+v0x*t,y0+v0y*t-0.5*g*t**2)

<IPython.core.display.Javascript object>

In [ ]: