Plotting EM diagrams from the flight model data?
-
Hi guys,
I’ve been reading the HFFM manual, and thought it would be useful/interesting to produce some EM diagrams of other aircraft.
I’ve written some matlab code to produce them, but although they seem approximately correct I think maybe they’re not quite right…
Could someone with knowledge of the flight models please give me a little guidance?
Here’s an example calculation for sustained turn (Ps=0) at sea level, F16-52, full AB, 100% fuel, drag index = 0, at mach 0.7
Density, rho = 0.0023769 slug/ft^3 Area, S = 300 ft^2 SpeedSound, a = 1116.45 ft/sec Weight, W = 19635 + 7162 = 26797 lb Gravity, G = 32.2 ft/sec^2 feet/sec, v = M a = 0.7 * 1116.45 = 781.515 Thrust, T = 31500 lbf = Drag = 0.5 * rho * S * v^2 * Cd :. Cd = T / ( 0.5 * rho * S * v^2 ) = 31500 / ( 0.5 * 0.0023769 * 300 * 781.515^2 ) = 0.1446551 Before looking up Cd on the Cd table, divide by its 0.66666 coefficient: Cd' = Cd / 0.66666 = 0.216985 Looking at the Cd table in the FM, at M=0.7, Cd'=0.216985 falls between the alpha breakpoints of 12.0 deg (Cd'=0.1965) and 13.0 deg (Cd'=0.2445) Interpolating linearly between those points gives: Alpha = 12.4268 deg At M=0.7 and Alpha=10.3541 the Cl table gives (between 0.871 and 0.942): Cl = 0.9013 Lift, L = 0.5 * rho * S * v^2 * Cl = 0.5 * 0.0023769 * 300 * 781.515^2 * 0.9013 = 196267 G load, n = L / W = 196267 / 26797 = 7.324 Turn rate, d = G * ( n^2 - 1 )^0.5 / v = 32.2 * ( 7.324^2 - 1 )^0.5 / 781.515 = 0.2989 radians/sec = 17.13 degree/sec ```(Actually, I've lost some rounding precision during the example - my code gives a more accurate 17.21 deg/sec) 17.21 deg/sec isn't a ridiculous figure, but looking at other EM diagrams from other sources it looks like the correct turn rate should probably be closer to 16 deg/sec That 1 deg/sec makes a difference, and gives my curves some peaks and dips that I don't think should be present. Could someone familiar with this please look at this example and tell me: A. Whether the technique is correct? B. Whether you come to the same result? Also, can anyone tell me about the Stores Drag limiter? It has these values: 0.900000 0.000100 1.000000 0.000283 Which I assume is two pairs of numbers: (0.9, 1.0) and (0.001, 0.000283) meaning when "XXX" = 0.9, Stores Drag = 1.0, and when "XXX" = 0.001 Stores Drag = 0.000283, with linear interpolation between. My problem is I don't know what "XXX" is! What variable are the two left hand numbers referring to? -
You can open dat files with this and it is able to create plots.
http://www.mediafire.com/download/z0bs9jg7217he83/f4doghouse.zipHave fun!
-
dont multiply CD by 0.66666 , this coefficient is just there because there is hack in the code that divide again by 0.666
so consider the coefficient is 1.0
-
You can open dat files with this and it is able to create plots.
http://www.mediafire.com/download/z0bs9jg7217he83/f4doghouse.zipHave fun!
Thanks! I’ll have a good look at this, it may help me with my charts.
dont multiply CD by 0.66666 , this coefficient is just there because there is hack in the code that divide again by 0.666
so consider the coefficient is 1.0
Ok, I’ve changed that. For the same example above, I now get 14.51 deg/sec, which seems a little low?
Are there any other multipliers in the code, one for thrust for example?Could you guys shed any light on the Stores Drag limiter for me?
-
btw rho ans Speed of sounds are of course variable
standard atmosphere is used
if you want to test your code, just look at HFFM EM charts that have been plotted from Dat files
so i advise you run your code on exact HFFM values
like
DI 38
GW 28330 -
btw rho ans Speed of sounds are of course variable
standard atmosphere is used
if you want to test your code, just look at HFFM EM charts that have been plotted from Dat files
so i advise you run your code on exact HFFM values
like
DI 38
GW 28330Thanks Mav, yep I’m using the atmosphere data from the HFFM manual.
I would very much like to compare mine exactly with the HFFM charts, but I can’t figure out the Stores_Drag coefficient limiter from the FM file - could you tell me what other variable it’s linked to?
For example in the block 52 the limiter is: 0.900000 0.000100 1.000000 0.000283
And I’m pretty sure that is two data points: (0.9, 1.0) and (0.001, 0.000283) where 1.0 and 0.000283 are the Stores_Drag coefficient, but what are 0.9 and 0.001 referring to? Mach? Cd? -
for OFM
CD += 0.9 * 0.0001 * DI if below mach 0.9
CD += 1.0 * 0.000283 * DI above mach 1.0
interpolate between
for AFM, values are stores in AFM file explicitely
-
Thanks very much Mav!
I can now compare against the HFFM manual. Unfortunately, although my charts are very similar, they aren’t the same
My charts on the left, HFFM manual on the right.
These show the maximum performance envelope in blue, and the sustained turn in green.
Whenever I plot sustained turn for an F16 I get that “bump” at around mach 0.6 which isn’t on the HFFM charts… I can’t understand it because the calculation at 0.6 is done precisely the same as everywhere else on the curve. I don’t get it for any other airframes besides the F16.
The righthand part (greater than mach 1.0) of my sustained turn curve is “stretched” along the x axis, easily noticeable by comparing the maximum speeds: mine are too high. I wonder if there’s some piece of code that I don’t know about that handles sub-sonic and super-sonic slightly differently?
My application of the alpha vs g limiters to the stall curve seems to be slightly too harsh…comparing the maximum instantaneous turn rate and the minimum turn radius mine are very slightly worse in both examples (my stall line near the peak is slightly more to the right, and slightly lower). So it looks like I have a problem there too.
Mav-JP, do you have any suggestions just from looking at these images?
You can open dat files with this and it is able to create plots.
http://www.mediafire.com/download/z0bs9jg7217he83/f4doghouse.zipHave fun!
I’ve had a look at this, it’s fantastic! So much potential. But I can’t see any way to plot Ps lines on a doghouse plot, which is a real pity - is there a way to do this?
And it doesn’t seem apply the aoa vs g limiters to the F16, which is a little unfortunate. -
Ah, oops, found the problem that was giving a too high top speed. The Ps=0 curve above mach 1.0 is perfect now.
The stall line is still applying the g limiter very slightly too harshly for some reason. My minimum radius is 1603 ft vs the correct 1546, and my max instantaneous turn is 21.2 deg/sec instead of the correct 21.4
And the Ps=0 line is fairly different below mach 1.0, using the following equations (includes DI):
Mach 0.7, sea level, F16bk52, max AB, DI = 38, 28330 lb Density, rho = 0.0023769 slug/ft^3 Area, S = 300 ft^2 SpeedSound, a = 1116.45 ft/sec Weight, W = 28330 lb Gravity, G = 32.2 ft/sec^2 feet/sec, v = M a = 0.7 * 1116.45 = 781.515 Thrust, T = 31500 lbf = Drag = 0.5 * rho * S * v^2 * Cd :. Cd = T / ( 0.5 * rho * S * v^2 ) = 31500 / ( 0.5 * 0.0023769 * 300 * 781.515^2 ) = 0.1446551 Cd = Cd0 + StoresDragCoeff * DI :. Cd0 = 0.1446551 - 0.9 * 0.0001 * 38 = 0.1412351 Looking at the Cd table in the FM, at M=0.7, Cd0=0.1412351 falls between the alpha breakpoints of 12.0 deg (Cd0=0.131) and 13.0 deg (Cd0=0.163) Interpolating linearly between those points gives: Alpha = 12.3198 deg At M=0.7 and Alpha=12.3198 the Cl table gives (between 0.871 and 0.942): Cl = 0.8937 Lift, L = 0.5 * rho * S * v^2 * Cl = 0.5 * 0.0023769 * 300 * 781.515^2 * 0.8937 = 194613 G load, n = L / W = 194613 / 28330 = 6.8795 Turn rate, d = G * ( n^2 - 1 )^0.5 / v = 32.2 * ( 6.8795^2 - 1 )^0.5 / 781.515 = 0.2800 radians/sec = 16.04 degree/secAgain, there is some slight rounding error there, my code gives 16.12. 16.12 is pretty good, pretty close, but the HFFM manual chart says it should actually be about 16.6 or 16.7….
What is wrong with my calculation above? Is the g value used in the code exactly 32.2?
This shows my curves in pink and green over the HFFM manual figure.
Also, I have a small query about weights in the HFFM manual. The EM diagrams are all plotted with 2 AIM9M + 2 AIM120 + ALQ131 attached, but according to the total weights in the manual, and the weights in the FM files, the payload weight varies? See below: the bottom line values should all be the same…?
What is the reason for this? Apologies if it’s a silly question, but I can’t find the answer.[b] bk50 bk52 bk32 bk30[/b] ----------------------------------------------------------------------- Dry weight from .dat file, D 19900 19500 18000 18700 Max. Fuel from .dat file, F 7162 7162 7162 7162 HFFM Manual EM weight, W 28670 28330 26609 27395 2 AIM9M + 2 AIM120 + ALQ131, = W - D - F [b] 1608 1668 1447 1533[/b] <-- why are these different? -
This is starting to become interesting …… :=)
-
@A.S:
This is starting to become interesting …… :=)
I fear that without a BMS FM expert like Mav-JP to review my posted calculation I can go no further
I don’t know what’s wrong with it. -
#define GRAVITY 32.177F
-
Ah, oops, found the problem that was giving a too high top speed. The Ps=0 curve above mach 1.0 is perfect now.
The stall line is still applying the g limiter very slightly too harshly for some reason. My minimum radius is 1603 ft vs the correct 1546, and my max instantaneous turn is 21.2 deg/sec instead of the correct 21.4
And the Ps=0 line is fairly different below mach 1.0, using the following equations (includes DI):
Mach 0.7, sea level, F16bk52, max AB, DI = 38, 28330 lb Density, rho = 0.0023769 slug/ft^3 Area, S = 300 ft^2 SpeedSound, a = 1116.45 ft/sec Weight, W = 28330 lb Gravity, G = 32.2 ft/sec^2 feet/sec, v = M a = 0.7 * 1116.45 = 781.515 Thrust, T = 31500 lbf = Drag = 0.5 * rho * S * v^2 * Cd :. Cd = T / ( 0.5 * rho * S * v^2 ) = 31500 / ( 0.5 * 0.0023769 * 300 * 781.515^2 ) = 0.1446551 Cd = Cd0 + StoresDragCoeff * DI :. Cd0 = 0.1446551 - 0.9 * 0.0001 * 38 = 0.1412351 Looking at the Cd table in the FM, at M=0.7, Cd0=0.1412351 falls between the alpha breakpoints of 12.0 deg (Cd0=0.131) and 13.0 deg (Cd0=0.163) Interpolating linearly between those points gives: Alpha = 12.3198 deg At M=0.7 and Alpha=12.3198 the Cl table gives (between 0.871 and 0.942): Cl = 0.8937 Lift, L = 0.5 * rho * S * v^2 * Cl = 0.5 * 0.0023769 * 300 * 781.515^2 * 0.8937 = 194613 G load, n = L / W = 194613 / 28330 = 6.8795 Turn rate, d = G * ( n^2 - 1 )^0.5 / v = 32.2 * ( 6.8795^2 - 1 )^0.5 / 781.515 = 0.2800 radians/sec = 16.04 degree/secAgain, there is some slight rounding error there, my code gives 16.12. 16.12 is pretty good, pretty close, but the HFFM manual chart says it should actually be about 16.6 or 16.7….
What is wrong with my calculation above? Is the g value used in the code exactly 32.2?
This shows my curves in pink and green over the HFFM manual figure.
http://oi58.tinypic.com/qs4net.jpgAlso, I have a small query about weights in the HFFM manual. The EM diagrams are all plotted with 2 AIM9M + 2 AIM120 + ALQ131 attached, but according to the total weights in the manual, and the weights in the FM files, the payload weight varies? See below: the bottom line values should all be the same…?
What is the reason for this? Apologies if it’s a silly question, but I can’t find the answer.[b] bk50 bk52 bk32 bk30[/b] ----------------------------------------------------------------------- Dry weight from .dat file, D 19900 19500 18000 18700 Max. Fuel from .dat file, F 7162 7162 7162 7162 HFFM Manual EM weight, W 28670 28330 26609 27395 2 AIM9M + 2 AIM120 + ALQ131, = W - D - F [b] 1608 1668 1447 1533[/b] <-- why are these different?why this assumption that EM charts are plotted with full fuel ?
there is a reason why “strange weights” and “strange DI” are used in HFFM manual…think about it
-
btw…careful to take correctly CAT1 limiter into account…
for a given AOA, you shall limit your G’s…
for instance if you find AOA = 22 from CD curve
then compute Gs from CL
then G= min (G, CAT1(22))
in that case G = 5.2 MAX
-
This post is deleted! -
#define GRAVITY 32.177F
Thanks!
why this assumption that EM charts are plotted with full fuel ?
there is a reason why “strange weights” and “strange DI” are used in HFFM manual…think about itBut all the EM charts in the HFFM manual say they’re plotted with 100% fuel?
btw…careful to take correctly CAT1 limiter into account…
Yes, I’m applying the limiter to both the stall line and to the Ps curves, although it doesn’t appear to be hit in the Ps curves.
In the stall line I need to find the g loading, n, for maximum lift for a velocity, which is a function of both alpha and n, because of the limiter.
so: n = fn( alpha, n )
I tried iteratively converging to find n, but it wouldn’t converge for all points, so now I use this method for each velocity:- Start at the minimum alpha in the limiter, 15 deg, find the lift and therefore the n value.
- Compare this n value to the limiter n value.
- Increase alpha slightly, and check again.
- Keep checking until I find an alpha where the lift n value is higher than the limiter n value, and then I use the alpha value preceding this one.
So I’m using the highest possible alpha that doesn’t exceed the g limiter.
It seems logical to me, and it works quite well but my stall line is slightly more conservative than the HFFM manual: my minimum radius is about 50 feet more, and my maximum turn rate is about 0.2 deg/sec less.
Is that the correct method for calculating the stall line with the g limiter? Is there a better method?
Thank you very much for your help so far, I really appreciate it. If you have an opportunity would you be able to look at this example calculation for a point on the Ps=0 curve and tell me if it is correct?
Mach 0.7, sea level, F16bk52, max AB, DI = 38, 28330 lb Density, rho = 0.0023769 slug/ft^3 Area, S = 300 ft^2 SpeedSound, a = 1116.45 ft/sec Weight, W = 28330 lb Gravity, G = 32.177 ft/sec^2 feet/sec, v = M a = 0.7 * 1116.45 = 781.515 Thrust, T = 31500 lbf = Drag = 0.5 * rho * S * v^2 * Cd :. Cd = T / ( 0.5 * rho * S * v^2 ) = 31500 / ( 0.5 * 0.0023769 * 300 * 781.515^2 ) = 0.1446551 Cd = Cd0 + StoresDragCoeff * DI :. Cd0 = 0.1446551 - 0.9 * 0.0001 * 38 = 0.1412351 Looking at the Cd table in the FM, at M=0.7, Cd0=0.1412351 falls between the alpha breakpoints of 12.0 deg (Cd0=0.131) and 13.0 deg (Cd0=0.163) Interpolating linearly between those points gives: Alpha = 12.3198 deg At M=0.7 and Alpha=12.3198 the Cl table gives (between 0.871 and 0.942): Cl = 0.8937 Lift, L = 0.5 * rho * S * v^2 * Cl = 0.5 * 0.0023769 * 300 * 781.515^2 * 0.8937 = 194613 G load, n = L / W = 194613 / 28330 = 6.8795 Turn rate, d = G * ( n^2 - 1 )^0.5 / v = 32.177 * ( 6.8795^2 - 1 )^0.5 / 781.515 = 0.2802 radians/sec = 16.06 degree/sec -
OOOOOOOOOOOk
first we need to compare apples and apples
so in order to test your routine and compare with HFFM manual you MUST use HFFM files from HFFM package (OFM) and NOT BMS / AFM files since calculations are not exactly the same
secondly, FROM HFFM files you should take LEF into consideration in CL (not CD) by
alpha = alpha + lefangle
lefangle = (float)min ((af->auxaeroData->lefMaxMach - af->mach) / 0.2F, 1.0F) * af->auxaeroData->lefMaxAngle * DTR;
BE AWARE that CL table in HFFM files have LEF substracted IIRC damned i am not so sure anymore , it’s been 10 years ALREADY LOL…anyway your difference might come from this…
anyway, FIRST start from HFFM files, not BMS files since computation has changed for OFM / AFM / lef stuff…
HFFM manual has been drawn with matlab routines similar to yours so you should get the same results
-
OOOOOOOOOOOk
first we need to compare apples and apples
so in order to test your routine and compare with HFFM manual you MUST use HFFM files from HFFM package (OFM) and NOT BMS / AFM files since calculations are not exactly the same
secondly, FROM HFFM files you should take LEF into consideration in CL (not CD) by
alpha = alpha + lefangle
lefangle = (float)min ((af->auxaeroData->lefMaxMach - af->mach) / 0.2F, 1.0F) * af->auxaeroData->lefMaxAngle * DTR;
BE AWARE that CL table in HFFM files have LEF substracted IIRC damned i am not so sure anymore , it’s been 10 years ALREADY LOL…anyway your difference might come from this…
anyway, FIRST start from HFFM files, not BMS files since computation has changed for OFM / AFM / lef stuff…
HFFM manual has been drawn with matlab routines similar to yours so you should get the same results
Thanks Mav, this information is gold!
Am I understanding this correctly: To compare to the HFFM manual I should use FM files from FF4+ or OF4.1+ only, and not any FM files from BMS4?
Is that correct?Also, what is DTR above? Turn rate?
I am going to have to put my work on this on hold for a couple of weeks, I have to finish my thesis by 8 Dec…
-
correct n compare with files from OF…
DTR = Degree to Radian
-
ooops sorry
F16 OFM has a AOA lef angle , limited to mach 1.0
LEFangle = (float)max (min(af->alpha,af->auxaeroData->lefMaxAngle)* DTR, 0.0f)
and
lefFactor =LEFangle / ( DTR * auxaeroData->lefMaxAngle);
AFM lef is more complex, but you dont need it at all
