Meteorology 341 Exercise 6: Trajectory Model Due October 15

Although air is a fluid, the motion of an air parcel is determined by the forces acting on it, just like the motion of a particle. Air parcels are subjected to a variety of forces, including gravity, the pressure gradient force, buoyancy, the Coriolis force, and frictional or turbulent drag. In the last exercise, we utilized the nearly exact balance (called hydrostatic balance) between the vertical component of the pressure gradient force and the buoyancy force to calculate the height of pressure levels in the atmosphere.

In this exercise, we will calculate the trajectory, or path, of a particle acted upon by gravity only. In later exercises, we will extend this trajectory model to include some of the other forces that act upon air parcels.

Consider the motion of a particle.The particle’s trajectory is described by the particle’s position as a function of time t. In our case, we will restrict our attention to the particle’s height z for simplicity. Then the trajectory is described by z(t).

How are trajectories related to forces? The force per unit mass is the acceleration a. By definition,

a = \frac{dw}{dt} ,  (1)

while

w = \frac{dz}{dt}  (2)

is the vertical velocity. This pair of ordinary differential equations relates trajectories to forces.

Consider the following three cases.

1. If a = 0, then w(t) = w(0) and z(t) = z(0) + w(0)t.

2. If a = a constant, then w(t) = w(0)+at and z(t) = z(0)+w(0)t+ \frac{1}{2}at^{2}.

3. If a = a(z, w), then z(t) can be obtained by numerical integration.

We will now describe a method to numerically integrate Eqs. (1) and (2). First consider Eq. (2). Formally integrate this over a time interval Δt to obtain

z(t + Δt) − z(t) = ∫_{t}^{t+Δt} w(t) dt.

Approximate the integral using the trapezoidal rule:

∫_{t}^{t+Δt} w(t) dt ≈ \frac{1}{2}[w(t) + w(t + Δt)]Δt.

This becomes more accurate as Δtbecomes smaller. An approximate solution of Eq. (2) is therefore

z(t + Δt) = z(t) + \frac{1}{2}[w(t) + w(t + Δt)]Δt.  (3)

Similarly, an approximate solution of Eq. (1) is

w(t + Δt) = w(t) + \frac{1}{2}[a(z(t), w(t)) + a(z(t + Δt), w(t + Δt))]Δt. (4)

Note that the acceleration may depend on z and w, and that the functional form of this dependence (that is, a(z, w)) will always be specified. Also note that z(t) and w(t) are considered to be known, while z(t+Δt) and w(t+Δt) are unknown. As a result, Eqs. (3) and (4) are a pair of linear algebraic equations and can be solved for the two unknowns. However, the algebra can often get quite complicated. With a simple approximation, we can convert Eqs. (3) and (4) into ones in which the unknowns depend explicitly on the knowns, rather than implicitly. We will replace z(t + Δt) and w(t + Δt) where they appear on the right-hand sides of Eqs. (3) and (4) with the following approximations:

z(t + Δt)^{∗} = z(t) + w(t)Δt

and

w(t + Δt)^{∗} = w(t) + a(z(t), w(t))Δt.

Then Eqs. (3) and (4) are replaced by

z(t + Δt) = z(t) + \frac{1}{2}[w(t) + w(t + Δt)^{∗}]Δt  (5)

and

w(t + Δt) = w(t) + \frac{1}{2}[a(z(t), w(t)) + a(z(t + Δt)^{∗}, w(t + Δt)^{∗})]Δt. (6)

Eqs. (5) and (6) can be written in a compact form as

z(t + Δt) = z(t) + Δz

and

w(t + Δt) = w(t) + Δw.

How do we use Eqs. (5) and (6) to calculate (numerical approximations to) z(t) and w(t) for 0 ≤ t ≤ T? Here is the algorithm:

1. At time t = 0, z(t) and w(t) are known.

2. Calculate z(t + Δt)^{∗} and w(t + Δt)^{∗}.

3. Calculate Δz and Δw.

4. Calculate z(t + Δt) = z(t) + Δz and w(t + Δt) = w(t) + Δw.

5. Replace t with t + Δt, z(t) with z(t + Δt), and w(t) with w(t + Δt).

6. If t < T, then go to step 2. Otherwise, stop.

A FORTRAN version of this algorithm:

t = 0.
z = z_0
w = w_0
do i = 1, imax
zstar = ...
wstar = ...
dz = ...
dw = ...
z = z + dz
w = w + dw
t = t + dt
end do

Assignment

1. Write a C or FORTRAN program that calculates a trajectory using the algorithm described above.

2. Use your program to calculate the trajectory of a particle acted upon by gravity only. Set a(z, w) = −g = −9.8 m s^{−2}. At t = 0, set z = 0 and w = 98 m s^{−1}. Set Δt = 0.2 s and T = 20 s. Plot the numerical and the exact solutions for z on one graph and for w on a seperate graph.