Naval Gun

And now a brief intermission from all that "space stuff". I thought I might tackle the challenge of modeling a ballistic trajectory of a 15" naval gun, as used on this ship, for example. Now there are a couple of factors that influence the trajectory of a projectile of a given velocity; gravity and drag.

Gravity I have already looked into in my previous posts. For this model, fired in a local area on the surface of the Earth, the variations in its acceleration would be negligible, so I established it constant at 9.807 m/s². At this point I will just go ahead and plot a trajectory without drag (air resistance) to see what I get.

Apart from the gravity, I established some initial data for the calculation. They would be shell velocity (initial, or muzzle velocity if you like), and angle (elevation) of the gun, so that it fires upwards like a mortar shot, the way it did, and loft the shell off towards its target in the distance. Here's how the program started off (I am not going to do that complete color formatting anymore, it is a pain. From now, comments are going to be pale blue, code in dark blue);

int main()
{
// variable declarations...
double grav_accel = 9.807;
double velocity = 785; // meters per second
double xpos = 0, ypos = 0; // The shell positions, each unit a meter.
double xvel, yvel; // The shell X and Y (2D) velocity components.
double angle = 0.40; // Angle in radians, approximately 23º.
// For the record (I am always mixing this up)....
// To convert radians to degrees, first divide them by PI, then multiply by 180.
// To convert degrees to radians, first divide them by 180, then multiply by PI.

// Code continues after the following comment....

Once that is established, we can set up a loop to move the shell. The thing to consider here is that during the upward travel of the shell, gravity is going to slow its velocity on the Y component (vertical). This will slowly alter the trajectory, as it will also slow its total velocity. And on the "way down" gravity will increase the projectile velocity (drag being nonexistent, to start with. Later we will see how drag affects this).

This is what needs to be calculated.

// Code continues here...

while(ypos >= 0) // While the shell is above the "ground", or "sea", as it were.
{
xvel = cos(angle) * velocity;
yvel = sin(angle) * velocity;

xpos += xvel;
ypos += (yvel - grav_accel);

velocity = sqrt(pow(xvel, 2) + pow((yvel - grav_accel), 2));
angle = atan2f((yvel - grav_accel), xvel);

usleep(1000000);
}
return 0;
}

At first glance it may appear that there is a bit much there. Why not just compute the vertical velocity component and leave it at that? There would be no change to the horizontal component of the velocity as it changed its trajectory, which is what really happens. That is why. As the shell travels upwards, a component of gravity is subtracting velocity from the shell's total velocity.

For example, if the code were left simply as...

xvel = cos(angle) * velocity;
yvel = sin(angle) * velocity;

xpos += xvel;
ypos += (yvel - grav_accel);

...and looped through as that, the horizontal velocity would always be the result of Cos(angle) x velocity (ie; Cos(22º) x 785 = 727.8 m/s). This would extend the range of the gun unrealistically, as there would be no degradation of its velocity caused by it "fighting its way uphill" at the beginning of its trajectory. The new velocity is calculated by invoking Pythagoras with the new components of X and Y velocity. Then the new angle of the resultant trajectory is also calculated from the new components (atan()), and that data is fed back into the loop. The program will run as it is, now.

However, that is not the end of the story. The shell travels through the atmosphere, and the atmosphere has a density that causes aerodynamic drag. The formula for drag is;

Drag = 0.5 x Velocity² x Air_Density x Drag_Coefficient x Frontal_Area

In SI units, that would give you the force of drag in Newtons (I believe. Most of my usage of lift and drag formula in the past have been in Imperial units, which provides the force in pounds. This is the first time I am applying it with SI units). So, now I needed to supply some more data to the program in order to calculate this. Velocity we already have, so we add to the variable declarations the following;

double frontal_area = 0.1148; // m², about the frontal area of a 15 inch diameter shell.
double Cd = 0.22; // An estimated, "probable" coefficient of drag.
double Rho = 1.225; // kg/m³, air density, according to SL ISA.
double drag;

All very well, but now the drag has to "work" against something. That "something" would be the momentum (inertia) of a shell of a certain mass traveling at a certain velocity. Here is where I get a bit "plauged" by doubt, but I do believe that I am correct in stating that the motion inertia of the shell is equivalent to its mass multiplied by its velocity. Let's call it "momentum" for the moment, to differentiate it from kinetic energy, which is a different formula (and provides preposterously large results for the purpose of this model). The unit of the product would be the Kg Meter per Second (Newton Second). I add these two more variables to the declarations in the code;

double mass = 870; // kg
double momentum;

And now to the calculations. This would go inside the loop, in the code above, just after the angle calculation and before the usleep() function call.

drag = 0.5 * Rho * pow(velocity, 2) * Cd * frontal_area;
momentum = (velocity * mass) - drag; // I believe the units are compatible here.
// At least, the results of the running program are quite acceptable.

velocity = momentum / mass;
// From the energy, devise the new velocity, and then loop back to the beginning,
// armed with the new trajectory angle and shell velocity.

It might not be perfect, but the output (which can be seen by adding some cout statements to the code) is actually quite convincing. There is a slight excess of range, but I believe it might be because the Cd is a tad low. After all, it was only estimated. The shell goes up almost 4 km, travels 29 km, taking 55 seconds to do so, and arrives at the target at a velocity of about 480 m/s, having decelerated from its original 785 m/s, and been at about 540 m/s at the top of the trajectory. These figures are not far off the performance of the real gun.

After note:
I had some fun with this one. Obviously, the code can be cleaned up, or neatly packaged away in a function, but the basic idea is there. A couple of members of a naval history forum provided me with some additional sources of information (data) on gun ranges and trajectories. I really was not far off, at all, except that it is closer to the real performance of the gun if the coefficient of drag is raised to 0.25 or 0.26, instead of 0.22.

'Till the next installment!