Plotting a Trajectory for a Ballistic Object: Complexities Revealed

[Zaggy]

OK, having been researching this of late, I've come to the conclusion it is far more complex than I first thought! So I guess without further ado, I'll jump into explaining my problem...

Aim:

Having established a set of initial conditions, I would like to be plot the trajectory of a ballistic object (in two dimensions, vertical and horizontal) at a repeating interval (say every 1/100th or 1/1000th of a second), taking into account gravity (obviously), accurate drag model & ballistic coefficent and the effects of an atmosphere (initially, the standard atmosphere)

Initial Conditions:

i) the 'world' / standard atmosphere
We can assume a standard world where g will be constant, though RHO (and speed of sound) may have a different value OR be changing during flight (ie as a projectile is shot upwards into thinner/colder air):

g (gravity) : -32.15 ft/s^2
h (altitude) : 0 ft
T (Temperature) : 59deg F (where T = 59 - 0.00356 * h)
p (pressure) : 2118.145 lbs/sqft (where p = 2116 * (((T + 459.7) / 518.6) ^ 5.256))
RHO (density) : 0.0023769 slugs/cuft (where RHO = p / (1718 * (T + 459.7)))
Speed of Sound : 1119.28 ft/s (assuming 80% humidity)​

ii) projectile properties
We are assuming a spin-stabilised projectile, with a circular cross-section (I'm trying to use a real-world set of numbers, based up high-power target shooting projectiles) that is fired horizontally:

v (velocity) : 3000 ft/s
W (weight) : 180grains
r (radius of cross-section) 0.154 inches (ie, a .308" cal projectile)
BC (Ballistic Coefficent) : 0.400 lbs/in^2 (G7 Drag Model - see excerpt)​

'G7' Drag Model excerpt (30 data points, instead of 80 - Drag Model assumes the Standard Atmosphere is present; so as I understand it, we don't even really need h, T, p and RHO; just a means to convert the current v (ft/s) into Mach; which is not difficult):

Mach CDrag
0.900 0.1464
0.925 0.1660
0.950 0.2054
0.975 0.2993
1.000 0.3803
1.025 0.4015
1.050 0.4043
1.075 0.4034
1.100 0.4014
1.150 0.3955
1.200 0.3884
1.300 0.3732
1.400 0.3580
1.500 0.3440
1.600 0.3315
1.700 0.3209
1.800 0.3117
1.900 0.3042
2.000 0.2980
2.100 0.2922
2.200 0.2864
2.300 0.2807
2.400 0.2752
2.500 0.2697
2.600 0.2643
2.700 0.2588
2.800 0.2533
2.900 0.2479
3.000 0.2424
​

Method:

Now, as above, the goal is to plot the movement of a project over time in 2D space. Plotting a trajectory with ONLY gravity is easy - plotting a trajectory with a CONSTANT retardation is easy (say 60% retardation per second, which is 'in the ballpark' for this sort of projectile) - BUT I have no idea about how to go about figuring out the retardation using the Ballistic Coefficent and the Drag model (obviously the gravity component is not a problem, just 'deceleration' or 'retardation' of the projectile)...

In otherwords, what do I need to in order to calculate velocity change from time= 0.000 to time = 0.001... eg: t = 0.000, v = 3000.0 --> t = 0.001, v = 3000.0 - x1 --> t = 0.002, v = (3000.0 - x1) - x2; etc...


I hope that makes sense to someone out there :)


Cheers


Zaggy

 

[rcgldr]

You need to implement some form of numerical integration. Below is an example of a predictor - corrector algorithm:

With the data you have, you can generate a function to calculate acceleration based on velocity and position:

Acceleration = F(velocity, position)

In the algorithm shown below, a_n, v_n, and p_n, are successive "guesses" that should converge quickly. a₁, v₁, and p₁ are Euler (since a1 is really a(i-1)), but I left the expressions similar to the ones used for trapezoidal (this can be used to simplify the code). The remaining steps are trapezoidal (average acceleration = 1/2 (previous acceleration + new acceleration)) and should converge quickly in 4 to 8 steps.

F(...) calculates the acceleration based on v and p, 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. 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).

You'll need to implement this algorithm for each axis.

initial condition:

time = 0
i = 0
v(0) = initial_velocity
p(0) = initial_position
a(0) = F(v(0), p(0))

then perform this loop until desired elasped time or impact

v₀ = v(i)
p₀ = p(i)

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

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

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

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

...

i += 1
time += Δt

v(i) = v_n
p(i) = p_n
a(i) = F(v(i), p(i))