More for the Naval Gun

As the title says, a bit more for the Naval Gun. I am going to do this in a couple of posts, to keep it in bite sized portions. Ballistics, especially of a spinning shell, is a complex, 3 dimensional arena. This model is limited to 2D for the moment. The day may come when I decide to tackle precession induced drift. For the moment, I just want to make this one, as it is, a bit more credible. One of the features of a naval shell, as fired from those magnificent warships of World War I and II, was that they were supersonic, a while before aircraft broke the legendary "sound barrier".

One of the notable traits of approaching "the speed of sound" is that the behavior of drag changes. There is a sharp increase in drag through the transonic region, and then a gradual fall off of this increased drag as the body accelerates beyond the speed of sound. This behavior is best represented, in a simulator model, by making the Coefficient of Drag a variable, the value of which is dependent upon the shell's velocity in relation to Mach 1. In my previous model, Cd was a fixed value. This was okay for starters, but it is NOT correct.

Now, things get a bit more intricate. The speed of sound is itself dependent on the atmospheric temperature. If you are a pilot, or even a dedicated "simmer", you will probably have an ARC-1 or E6B Flight Computer. There is a little scale on the computer side, where you match up a temperature to an M symbol, and read off the speed of an equivalent Mach number from the inner scale to knots on the outer scale. As a rough guide, a velocity of 340.278 m/s is equivalent to Mach 1 at 15º C (International Standard Atmosphere at Sea Level). And every drop of 1º C incurs a reduction of Mach 1 equivalent velocity by 0.61 m/s, approximately. So, if the temperature is 12º C, then Mach 1 is equivalent to;

340.278 - ((15 - 12) * 0.61)

To get the Mach number from this, divide the actual velocity of the shell (in m/s) by the corresponding Mach 1, in m/s, obtained for the given temperature.

Mach_Number_of_Shell = Velocity_of_Shell_ms / Velocity_of_Mach_1_ms

And that, basically, is what I am going to do with my model as a required step in determining a variable Coefficient of Drag. The temperature drops at an environmental lapse rate of 1º C per every 152.39 meters gained.

Declarations first, with all the others;

double ISA_temp_C = 15;
double temp_meters_degree_C = 152.39; // Meters of altitude per every -1º C
double ref_Mach_1 = 340.278; // Meters per second at 15º C
double Mach_fall_off = 0.61; // Mach meters per second less per every -1º C
double velocity_Mach;

Now, I decided to make a small function that I would call from inside the update loop. Here it is;

double Mach_Calculator(double altitude, double SL_temp, double alt_temp_lapse, double ref_M1, double Mach_lapse, double vel)
{
double temp_at_alt = SL_temp - (altitude / alt_temp_lapse);
double Mach_1_at_alt = ref_M1 - ((SL_temp - temp_at_alt) * Mach_lapse);
return vel / Mach_1_at_alt;
}

And that was it, in the simplest terms. I now have the shell velocity expressed in Mach number. I call the function with these parameters;

velocity_Mach =
Mach_Calculator(ypos, ISA_temp_C, temp_meters_degree_C, ref_Mach_1, Mach_fall_off, velocity);

Onwards. Now let's go get that variable Coefficient of Drag. There are (very complicated) ways of computing Cd. I opted for an easy way out, at this stage. Check out this page. Down on the Doppler Radar-measurements section, there is a model of a bullet's coefficient of drag at different Mach velocities. I mapped something like that into two arrays, declared globally, in my program, like this;

const int MACH_CD_ELEMENTS = 9;

double ref_Mach_table[MACH_CD_ELEMENTS] =
{-0.10, 0.00, 0.40, 0.60, 0.87, 1.02, 1.15, 2.40, 3.40};

double ref_Cd_table[MACH_CD_ELEMENTS] =
{0.000, 0.231, 0.229, 0.172, 0.139, 0.333, 0.341, 0.268, 0.255};

Simply matched up some Mach numbers (first array) to some Cd values (second array. Now it will just be a case of using the computed Mach velocity of the shell to find a corresponding Cd. So, it is all very well to say, looking at the two arrays, that if the shell is traveling at Mach 1.15, then the Cd is 0.341 (both being the seventh element in their respective arrays), but what if the shell is doing Mach 1.95? The Cd value is somewhere between 0.341 and 0.268. We need to interpolate. Here's the function that does it...

double Cd_Calculator(double Mach)
{
int index_count;

for(index_count = 0; index_count < MACH_CD_ELEMENTS; index_count++)
{
if(ref_Mach_table[index_count] > Mach)
break;
}

double inter_factor =
(Mach - ref_Mach_table[index_count - 1]) / (ref_Mach_table[index_count] - ref_Mach_table[index_count - 1]);

double output =
((ref_Cd_table[index_count] - ref_Cd_table[index_count - 1]) * inter_factor) + ref_Cd_table[index_count - 1];

return output;
}

So, what's it doing? That's the subject of the next post, where I will also interpolate to get a variable atmospheric density at different altitudes, adding a bit more polish to the model. For the moment, this function gets called AFTER determining the Mach number and BEFORE calculating the drag (in the update loop), like this...

Cd = Cd_Calculator(velocity_Mach);

That's already improved the model. It now has a variable Coefficient of Drag, dependent on the Mach number, which is in turn dependent on the shell velocity at different air temperatures.