Velocity Dependent Acceleration within Drag Force

[palnm]

For fun, I wanted to crate a program that dealt with 2d projectile motion w/ wind resistance. Specifically, where the object started from rest at a given height. I looked up several equations that dealt with wind resistance. Specifically, drag force. Drag force was defined as Fd = CdPAV^2/2. So, I thought acceleration would be a = CdPAV^2/2*m. Dealing with just vertical acceleration the equation would read a = g - Fd/m. I thought that I could double integrate both the vertical and horizontal accelerations and find, parametrically, a general position v. position equation, for that object. The problem that was immediately apparent was that Fd was dependent on velocity. Since this is acceleration equation, how do I integrate velocity FOR velocity. I figure I am doing something wrong. Could someone explain what?

 

[rcgldr]

Note this post only describes the math for the y-axis. (The x-axis movement would not involve g).

This can be integrated directly to determine velocity versus time, but apparently not position versus time. Wiki article with the answers, but not the derivation:

http://en.wikipedia.org/wiki/Drag_(physics)#Velocity_of_a_falling_object

http://en.wikipedia.org/wiki/Terminal_velocity#Derivation_for_terminal_velocity

derivation in this pdf file (see section 3, falling object with Newton drag):

http://math.kennesaw.edu/~sellerme/sfehtml/classes/math2202/fallingbodies.pdf

You can use numerical integration using the equation for acceleration versus velocity, to calculate velocity and position over a series of small time steps, or use the equation linked to above to calculate positions using the equation for velocity versus time. Example using acceleration versus velocity:

Acceleration = F(v) = g - ( Cd ρ A v²) / (2 m)

For each time step, you calculate

new_velocity = old_velocity + time_step x average_acceleration

One way to do this is to use trapezoidal method

average_accleration = 1/2 (old_acceleration + new_acceleration)

new_velocity = old_velocity + time_step x 1/2 (old_acceleration + new_acceleration)

where old_acceleration is the acceleration at the start of a time step and new_acceleration is the acceleration at the end of a time step. This requires you be able to calculate new_acceleration directly or to be able to predict it via an iterative method (see below)

You can then calculate position once you calculate the new velocity, using the same method:

new_position = old_position + time_step x 1/2 x (old_velocity + new_velocity)

An iterative corrector method will improve the results. In the algorithm shown below, a_n, v_n, and p_n, are successive "guesses" that should converge quickly. F(...) calculates the acceleration based on v_n(t), and Δt is the elapsed time per step. You may want to do 6 to 8 interations instead of the 4 shown in this example. The first step is essentially Euler (since a₁(t) is set = F(v₀(t)) (= a(t-1)), the remaining steps are trapezoidal. Even though each step of this algorithm will converge to a specific set of values, the algorithm is based on trapezoidal rule, a linear approximation, so you need to use small time steps (Δt).

F(v) = g - ( Cd ρ A v²) / (2 m)

v₀(t) = v(t-1)
p₀(t) = p(t-1)

a₁(t) = F(v₀(t))
v₁(t) = v(t-1) + 1/2 (a(t-1) + a₁(t)) Δt
p₁(t) = p(t-1) + 1/2 (v(t-1) + v₁(t)) Δt

a₂ = F(v₁(t))
v₂(t) = v(t-1) + 1/2 (a(t-1) + a₂(t)) Δt
p₂(t) = p(t-1) + 1/2 (v(t-1) + v₂(t)) Δt

a₃ = F(v₂(t))
v₃(t) = v(t-1) + 1/2 (a(t-1) + a₃(t)) Δt
p₃(t) = p(t-1) + 1/2 (v(t-1) + v₃(t)) Δt

a₄ = F(v₃(t))
v₄(t) = v(t-1) + 1/2 (a(t-1) + a₄(t)) Δt
p₄(t) = p(t-1) + 1/2 (v(t-1) + v₄(t)) Δt

...

v(t) = v_n(t)
p(t) = p_n(t)
a(t) = F(p_n(t))

time += Δt
t += 1

This is a predictor-corrector type algorithm:

http://en.wikipedia.org/wiki/Predictor-corrector_method#Euler_trapezoidal_example

Since you have an equation for velocity versus time v = V(t), you can use it directly and just use trapezoidal rule (time steps still need to be small):

t = 0
time = 0
v(t) = v(time)
p(t) = p(time)

loop:

t += 1
time += Δt
v(t) = V(time)
p(t) = p(t-1) + 1/2 (v(t-1) + v(t)) Δt

 

[cepheid]

Welcome to PF,

I agree that it looks nasty: right now you've got what's called an integral equation. The velocity depends on its own integral in some way. My advice would be NOT to take that step of integrating both sides of the equation. Instead, leave it as it is. Then you have a differential equation -- the velocity is related to its derivative. That's something that I, at least, have a better idea of how one might go about solving. For instance, in the x-direction where drag is the only force acting, you have:

F_x = -\frac{1}{2}\rho A C_d v_x^2

ma_x = -\frac{1}{2}\rho A C_d v_x^2

You can write this as

m\frac{dv_x}{dt} = -\frac{1}{2}\rho A C_d v_x^2

This differential equation is separable.

 

[rcgldr]

v(t) = v_n(t)
p(t) = p_n(t)
a(t) = F(p_n(t))

time += Δt
t += 1

a(t) is a function of velocity, so that third line should be:

a(t) = F(v_n(t))