Free Fall with Drag and Varying Air Density

Question

I am trying to plot a graph in Excel describing a free fall of an iron sphere dropped from an altitude of 12 km, while taking into account the different levels of air density affecting the drag.

The function looks like this: $$v\left(t\right)=\frac{m\cdot g}{b}\cdot \left(1-e^{\frac{-b\cdot t}{m}}\right)$$ where $b$ is the the constant of air resistance (that considers the area and shape of the object along with the density of the medium). The problem is that throughout the fall, the constant isn't constant, as air density increases with lower altitude.

Therefore the function requires plugging in both $t$ and $b$ – which are mutually dependent. That is to say, I simply don't know what $b$ corresponds to a given time (other than the initial $t_0$ at 12 km).

Even to calculate the terminal velocity $v_T=\frac{m\cdot g}{b}$ seems to be a problem, as I can't simply plug in any one $b$.

Let's say I am okay with simply taking a few discrete table air density values. What am I missing here in order to be able to "compose" the final graph?

Answer 1

You need to go a step back and understand where the formula $$v(t)=\frac{mg}{b}\left(1-e^{-bt/m}\right) \tag{1}$$ comes from. Actually this formula (1) for $v(t)$ is the solution of the differential equation $$\frac{dv(t)}{dt}=g-\frac{b}{m}v(t), \quad\text{ with starting condition }v(0)=0 \tag{2}$$

When the drag value $b$ is constant, then this differential equation is simple enough to solve and find the solution (1).

But, as you already noticed, $b$ is not constant. Instead it depends on the height $h$, which in turn depends on the time $t$. A simple approximation is (assuming constant temperature across the whole atmosphere) $$b(h)=b_0 e^{-h/H} \tag{3}$$ where $b_0$ is the drag value at ground level ($h=0$) and $H$ ($\approx 8$ km) is the so-called scale height of the earth's atmosphere.

So now you need to use this $h$-dependent $b$ in differential equation (2). Furthermore, the differential equation for $v(t)$ is not enough. You also need a differential equation for $h(t)$. Finally you arrive at these two coupled differential equations. $$\begin{align} \frac{dv(t)}{dt}&=g-\frac{b(h(t))}{m}v(t), &\quad\text{ with starting condition }v(0)=0 \\ \frac{dh(t)}{dt}&=-v(t), &\quad\text{ with starting condition }h(0)=h_0 \end{align} \tag{4}$$ Unfortunately the differential equations (4) are too difficult to solve and to write down $v(t)$ and $h(t)$ as analytical functions similar as we could do in (1).

But you can still find an approximative solution by numerical methods, for example by Euler's method. Applying this method to the differential equations (4) you get the following iteration formulas. $$\begin{align} v(t+\Delta t)&=v(t)+\left(g-\frac{b(h(t))}{m}v(t)\right)\Delta t \\ h(t+\Delta t)&=h(t)-v(t)\ \Delta t \end{align} \tag{5}$$ You apply the formulas (5) over and over again:

• Start with $h(0)=h_0$ and $v(0)=0$ at time $t=0$
• From these you calculate $h(\Delta t)$ and $v(\Delta t)$ at time $t=\Delta t$
• From these you calculate $h(2\Delta t)$ and $v(2\Delta t)$ at time $t=2\Delta t$
• ...

In order to get a good approximation you should use a small value for $\Delta t$. I guess $\Delta t=0.1$ sec should be small enough.

Doing this calculation in Excel is certainly possible. But in my opinion it would be preferable to use a real programming language for these calculations.

Answer 2

If b is not constant, viz., $b\equiv b(h(t))$, then you will not get solution for $v(t)$ as given above. While solving the differential equation, it assumes $b=$ const. You have to solve differential equation, with $b=b(h(t))$. But, if you make assumptions, that density varying slowly with height, then you can use the Euler method to write velocity at a particular time. If the initial velocity is $v(t=0)=v_0$ at height $h_0$, then at small time interval $\Delta t$, the velocity would be

$$v(\Delta t) = \frac{mg}{b(h_{\Delta t})}.\Bigg(1-\exp\bigg(\frac{-b(h_{\Delta t})\Delta t} {m}\bigg)\Bigg)$$

Where you have to plug $h_{\Delta t} = h_0-v_0\Delta t$. At time $2\Delta t$, $$v(2\Delta t) = \frac{mg}{b(h_{2\Delta t})}.\Bigg(1-\exp\bigg(\frac{-b(h_{2\Delta t})2\Delta t} {m}\bigg)\Bigg)$$ this time you have to plug $h_{2\Delta t} = h_{\Delta t}-v(\Delta t)\Delta t$. And so on $\dots$

If you know the exact differential form the of drag equation, then it is easier to compute more precisely. Then we can use RK4 method to integrate the differential equation.

For calculating terminal velocity you can put on condition $$v(n\Delta t) \approx v((n-1)\Delta t)$$ Computationaly, it is if v(ndt)-v((n-1)dt) < 10^-100

In the above discussion, $b=6\pi\eta r$ for low Reynold number object in the fluid. The slightly more general expression for $b$ which takes the density of fluid into account is

$$b=\frac{1}{2}\rho vC_d A$$

where A is area and $C_d$ is drag coefficient. How $\rho$ change with height for low altitude is given by $$\rho = \rho_0 e^{-h/H_n}$$

Use all aforementioned expression to the plot.