Problem with drag in Ballistic Calculator

I have been working on this for a while and I know that it works without drag, a_x = 0 and a_y = gravity. And the model matches actual trajectories with correct coefficients when the AoS(angle of shot) = 0, but if the AoS is changed to any angle above or below 0 the projectile goes wild (falls faster than gravity or against gravity). I believe it is a problem with the coefficient of drag but I cannot find any solution or figure out what is going on. The problem is visible if AoS is changed to i.e. 10 or -10, graph one will show that the projectile is going wild. The additional graphs 2-10 at the bottom are for diagnostics but I cannot find anything out of the ordinary for them. Any help would be appreciated.

%Ballistic_Calculator

clc

clf

clear all

%%%vars%%%

cal = 5.56; %[mm] caliber

weight = 3.56; %[g] weight of bullet

V_initial = 980; %[m/s] initial velocity

BC = .505; %[lb/in^2] Balistic Coefficient standard

AoS = 0; %[deg] Angle of Shot

SIA = .0892; %[deg] Angle when Sighted In, sets zero distance ~(.0892)

AoB = AoS+SIA; %[deg] Angle of Bullet

Air_Density = 1.225; %[kg/m^3]

HOB = 2.56; %[inches] height over bore

gravity = -9.81; %[m/s^2] gravity

%%%Derived vars%%%

diameter = cal/1000; %[m] diameter of bullet

radius = diameter/2; %[m] radius of bullet

CS_area = pi*radius^2; %[m^2] cross sectional area of bullet

mass = weight/1000; %[kg] mass of bullet

BCa = BC*703.07; %[kg/m^2] Ballistic Coefficient Adjusted for metric

CD = mass/(BCa*CS_area); %[unitless] Coefficient of drag

HOBa = HOB*.0254; %[m] height over bore adjusted

%%%simulation vars%%%

T_max = .5; %[s] max time

dt = .001; %[s] time step time

t = [0:dt:T_max];

%%zero arrays%%

a_x = zeros(1,T_max*dt);

a_y = zeros(1,T_max*dt);

v_x = zeros(1,T_max*dt);

v_y = zeros(1,T_max*dt);

p_x = zeros(1,T_max*dt);

p_y = zeros(1,T_max*dt);

v = zeros(1,T_max*dt);

FD_x = zeros(1,T_max*dt);

FD_y = zeros(1,T_max*dt);

FD = zeros(1,T_max*dt);

a = zeros(1,T_max*dt);

%%initial positions%%

a_x(1) = 0;

a_y(1) = gravity; %[m/s^2]

v_x(1) = V_initial*cosd(AoB);

v_y(1) = V_initial*sind(AoB);

p_x(1) = 0;

p_y(1) = -HOBa;

v(1) = V_initial;

FD_x(1) = 0;

FD_y(1) = 0;

FD(1) = 0;

a(1) = 0;

%%Numerical Integration%%

for k = 1:length(t)

a_x(k+1) = FD_x(k)/mass; %%this seems to be the cause of the problem(s)%%

v_x(k+1) = v_x(k)+.5*dt*(a_x(k)+a_x(k+1));

p_x(k+1) = p_x(k)+.5*dt*(v_x(k)+v_x(k+1));

FD_x(k+1) = -.5*Air_Density*v_x(k)^2*CD*CS_area;

a_y(k+1) = a_y(1)+FD_y(k)/mass;

v_y(k+1) = v_y(k)+.5*dt*(a_y(k)+a_y(k+1));

p_y(k+1) = p_y(k)+.5*dt*(v_y(k)+v_y(k+1));

FD_y(k+1) = -.5*Air_Density*v_y(k)*abs(v_y(k))*CD*CS_area;

v(k+1) = sqrt(v_x(k)^2+v_y(k)^2);

a(k+1) = abs(sqrt(a_x(k)^2+a_y(k)^2)/gravity);

end

%%trimming arrays%%

p_x = p_x(1:end-1);

p_y = p_y(1:end-1);

v_x = v_x(1:end-1);

v_y = v_y(1:end-1);

a_x = a_x(1:end-1);

a_y = a_y(1:end-1);

FD_x = FD_x(1:end-1);

FD_y = FD_y(1:end-1);

v = v(1:end-1);

a = a(1:end-1);

%%%plots%%

figure(1)

subplot(5,2,1);

plot(p_x,p_y)

title('Trajectory')

hold on

x = 0:max(p_x);

LOS = (sind(AoS)/cosd(AoS))*x;

plot(x,LOS)

subplot(5,2,2);

plot(p_x,v)

title('Velocity')

hold on

plot([0,max(p_x)],[823,823]) %Effective Velocity

plot([0,max(p_x)],[343,343]) %Speed of Sound

subplot(5,2,3);

plot(t,p_x);

title('Distance over Time')

subplot(5,2,4);

plot(p_x,a);

title('Acceleration, Gs')

subplot(5,2,5);

plot(p_x,a_x)

title('Acceleration, x-axis')

subplot(5,2,6);

plot(p_x,a_y)

title('Accelaration, y-axis')

subplot(5,2,7);

plot(p_x,v_x)

title('Velocity, x-axis')

subplot(5,2,8);

plot(p_x,v_y)

title('Velocity, y-axis')

subplot(5,2,9);

plot(p_x,FD_x)

title('Drag Force, x-axis')

subplot(5,2,10);

plot(p_x,FD_y)

title('Drag Force, y-axis')

Answer 1

[David Goodmanson]

Hi John,

I wonder a bit about the units and dimensions, but let's say they are good. There is a problem with FD, as you surmise. Leaving aside the k --> k+1 updating aspect for the moment, you have

FD_x = -.5*Air_Density* v_x^2        *CD*CS_area;
FD_y = -.5*Air_Density* v_y*abs(v_y)* CD*CS_area;

These should be

absv = sqrt(v_x^2+v_y^2)
FD_x = -.5*Air_Density*  v_x*absv  *CD*CS_area;
FD_y = -.5*Air_Density*  v_y*absv  *CD*CS_area;

The drag force is proportional to v^2, and directed along -v. Its x component is

   -(v_x/|v|)*|v|^2 = -v_x*|v|.

Same for y. Note that because of the |v| factor, velocity in the y direction contributes to FD_x, and vice versa.

Also, the only term you have involving g is a(k+1), but it is not used to update the calculation (and its meaning is not clear). For a_y you need a second term, which is simply

   a_y = FD_y/mass -g.

It's good to work out how the time iteration goes, but this can all be done more accurately using a diffeq solver like ode45.

[John Goodling]

This solution worked flawlessly and is much appreciated. In the line a_y(k+1) = a_y(1)+FD_y(k)/mass, the a_y(1) is the initial condition and a_y(1) is "gravity" or = -9.81. A bit convoluted but I am the only one who has to use this code. ODE45 could have been used but I like to see the process and the difference in accuracy is negligible for me. This is just an exercise to sharpen and maintain my coding skills in a meaningful way.