Atmospheric drag

Part II: Numerical Methods

This page originally compiled Mar 2007 by TFR
One might think that there is a mathematical way of modeling anything, including extremely complex subjects like the motion of objects through a non-uniform atmosphere. Actually, even the most basic computations we are concerned with here, the combined forces of gravity and air drag on moving objects, are not easily modeled, even with calculus! The reason is that our simple technique of splitting motion into its horizontal and vertical components turns out to be not so simple after all. Atmospheric drag, unlike gravity, cannot be modeled in one dimension; the components are not independent, and changing one also changes the other.

For these problems, we turn to a technique used by engineers and scientists: numerical modeling. The basic idea is that if we calculate something, like the motion of an object, in small enough slices, we can approximate the action of the forces over that small time slice reasonably well. To approximate the actions over large time scales, we simply calculate many small successive time slices. Before the easy availability of computers, this meant many tedious hours with pen, paper, and logarithm tables. Such computations can now be set up relatively quickly in a spreadsheet or a programming language, and the computer never tires of repeating them.

Let us now return to our problem of the motion of a falling body in the atmosphere. This relatively simple problem, for which we have already found an exact mathematical ("analytical") solution, makes a good test case against which we can compare our numerical results.


First, we will compute a table of data using the analytical equations from "Atmosphere Part I",
    v = -vt tanh(gt / vt),   for velocity,
    a = -g(1 - v2 / vt2),   for acceleration, and
    s = (-vt2 / g) ln[cosh(gt / vt)]   for vertical position
For this exercise, we will use a baseball as our object, where:
    m = 0.142 kg
    diameter = 7.3 cm,   r = 3.65 cm = .0365 m,   A = πr2 = 0.00419 m2
    Cd = 0.3, a guess based on the shape of the object
Our physical constants will be set to:
    ρ = 1.2 kg/m3, an approximate average value of atmospheric density at low airspeed and altitude
    g = 9.81 m/s2, an approximate average value of gravitational acceleration at low altitude
So our terminal velocity calculates as:
    vt = √[ 2mg ] / [ ρACd ] = 42.98 m/s

As an example, let us calculate the first row of our table, corresponding to the first second of elapsed time:
    t = 1 (second)
    v = -vt tanh(gt / vt)   = -42.98 × tanh( (9.81 × 1) / 42.98)   = -9.643 (meters/second)
    a = -g(1 - v2 / vt2)   = -9.81 × (1 - 9.6432/42.982)   = -9.316 (meters/second2)
    s = (-vt2 / g) ln[cosh(gt / vt)]   = (-42.982/9.81) × ln[ cosh( (9.81 × 1) / 42.98) ]   = -4.863 (meters)

These results become the first row of our data table, and we compute the rest of the data, plugging the three above equations into our spreadsheet or programming language and repeating for as many seconds as we would like to continue:

 t	 v		 a		 s
 1	-9.643099031	-9.31610863	-4.862987976
 2	-18.36176219	-8.019283231	-18.97411723
 3	-25.55502799	-6.34142647	-41.07291964
 4	-31.05478578	-4.687815007	-69.5159553
 5	-35.01995607	-3.296275625	-102.669374
 6	-37.75930057	-2.237381712	-139.1471942
 7	-39.59644621	-1.482577434	-177.8878726
 8	-40.80413532	-0.966859156	-218.1310462
 9	-41.58763293	-0.623996898	-259.355427
10	-42.0915915	-0.400016162	-301.2136497
11	-42.41395965	-0.255327153	-343.4784451
12	-42.61944115	-0.162524481	-386.002854
13	-42.7501224	-0.103270914	-428.6925572
14	-42.83311324	-0.065546972	-471.4873081
15	-42.88576962	-0.041573791	-514.3487405

Some things to note: As our baseball falls, it approaches terminal velocity. We can see, in the second column, velocity approaching that value, 42.98 m/s. In the third column, acceleration approaches zero as drag increases to match gravity. In the first and fourth columns, we note that the baseball has fallen 514 meters in 15 seconds and is still not quite up to terminal velocity.
We note that the computer, typically, has given us data to absurd decimal places; our input data does not justify this level of accuracy, but that is how computers work.
This data will be our "gold standard". As it was calculated analytically, we assume it to be absolutely accurate, and we will want our numerical simlation to match it as closely as possible.


Next, we attempt to numerically simulate the above data using a much simplified set of equations. We will base our numerical simulation on the acceleration equation only, and approximate velocity and position by summation of the results.

This time we will compute the first three rows of the table, as an example of how to begin this sort of simulation. We start at t=1, v=0:

    t = 1 (second)
    a = -g(1 - v2 / vt2)   = -9.81 × (1 - 02/42.982)   = -9.81 (meters/second2)
    vinit = 0,   velocity at the start of this second
    vfinal = vinit + a   = -9.81 (meters/second),   velocity at the end of this second
    vavg = (vinit + vfinal)/2   = -4.905 (meters/second),   average velocity for this second
    s = vavg   = -4.905 (meters),   vertical position at the end of this second

    t = 2
    a = -g(1 - v2 / vt2)   = -9.81 × (1 - -9.812/42.982)   = -9.299,    where v is taken from vfinal of t=1
    vinit = -9.81,   taken from vfinal of t=1
    vfinal = vinit + a   = -19.109
    vavg = (vinit + vfinal)/2   = -14.459
    s = [s of t=1] + vavg   = -19.364

    t = 3
    a = -g(1 - v2 / vt2)   = -9.81 × (1 - -19.1092/42.982)   = -7.871,    where v is taken from vfinal of t=2
    vinit = -19.109,   taken from vfinal of t=2
    vfinal = vinit + a   = -26.980
    vavg = (vinit + vfinal)/2   = -23.044
    s = [s of t=2] + vavg   = -42.409

The above examples become the first three rows of our table of numerical data, the rest of the rows being simply repetition of the same technique used for t=2 and t=3 above:

t	a		v init		v final		v avg		s
 1	-9.81		0		-9.81		-4.905		-4.905
 2	-9.29886432	-9.81		-19.10886432	-14.45943216	-19.36443216
 3	-7.870597562	-19.10886432	-26.97946188	-23.0441631	-42.40859526
 4	-5.943974181	-26.97946188	-32.92343606	-29.95144897	-72.36004423
 5	-4.052837446	-32.92343606	-36.97627351	-34.94985479	-107.309899
 6	-2.548196971	-36.97627351	-39.52447048	-38.25037199	-145.560271
 7	-1.512823966	-39.52447048	-41.03729445	-40.28088246	-185.8411535
 8	-0.865509144	-41.03729445	-41.90280359	-41.47004902	-227.3112025
 9	-0.484237604	-41.90280359	-42.38704119	-42.14492239	-269.4561249
10	-0.267451246	-42.38704119	-42.65449244	-42.52076682	-311.9768917
11	-0.14664931	-42.65449244	-42.80114175	-42.72781709	-354.7047088
12	-0.080088453	-42.80114175	-42.8812302	-42.84118598	-397.5458948
13	-0.043641639	-42.8812302	-42.92487184	-42.90305102	-440.4489458
14	-0.023752435	-42.92487184	-42.94862427	-42.93674806	-483.3856938
15	-0.012919018	-42.94862427	-42.96154329	-42.95508378	-526.3407776

We can compare our numerical results with the analytical "gold standard" to check how well it worked; column "v final" with column "v" above, and columns "a" and "s". We note differences, namely most of the data is somewhat higher in magnitude in our numerical simulation. We can attribute this to the coarseness of our simulation. Had we chosen a smaller time interval, say 1/4 second, summing the data 4 times for each full second, the data would have better approximated the analytical. We must take care in doing so that we account for it in the summation of final velocity, dividing acceleration by 4 in this case, and in summation of position, dividing average velocity by 4 as well. Here is the data for t=1, calculated at 1/4 second intervals, for comparison:

t	a		v init		v final		v avg		s
0.25	-9.81		0		-2.4525		-1.22625	-0.3065625
0.50	-9.77805402	-2.4525		-4.897013505	-3.674756752	-1.225251688
0.75	-9.682631866	-4.897013505	-7.317671471	-6.107342488	-2.75208731
1.00	-9.525590565	-7.317671471	-9.699069113	-8.508370292	-4.879179883

The output data at 1 second is better, but still not exact. We can go as far as we like with this, to the limit of the ability of our computer program to repeat the calculations. Here is a comparison of our data at t=1 second, calculated with diminishing time steps:

step	a		v final		s
1	-9.81		-9.81		-4.905
0.5	-9.68221608	-9.74610804	-4.88902701
0.2	-9.487604302	-9.688562534	-4.8765645
0.1	-9.405588755	-9.666504851	-4.870454613
0.05	-9.361735543	-9.654967141	-4.866904414
0.02	-9.334566154	-9.647885502	-4.864600193
0.01	-9.325371339	-9.645498769	-4.863801831

analyt	-9.31610863	-9.643099031	-4.862987976

Note that the case of time step = 0.01 s is the summation of 100 calculations! It is obvious that the smaller the time interval and the greater the number of summations for each second of elapsed time, the closer to the analytical standard we come.

The point is that we can actually simulate the motion of our object with this relatively simple method to an arbitrary accuracy, and certainly good enough to estimate the fall of a baseball. What this means is that we can also simulate the motion of a rocket through the atmosphere well enough to estimate its actual flight path, using similar simple techniques.


Home


Last Revised: Mar 2007