Weird response in take-off run

[seanmcleod]

I'm planning on trying to calculate $V_{MCG}$ for a given aircraft by running a set of simulated $V_{MCG}$ test runs using JSBSim.

So the very initial code I wrote, see below, simply loads the 737 model from the JSBSim repo, sets the engines running, performs a ground trim and then advances the throttles to max and then graphs the airspeed as a function of time for 120s.

import jsbsim
import matplotlib.pyplot as plt

def plot(x_data, y_data, title, xlabel, ylabel):
    plt.plot(x_data, y_data)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.title(title)
    plt.show()

# The path supplied to FGFDMExec is the location of the folders "aircraft", "engines" and "systems"
fdm = jsbsim.FGFDMExec('..\\') 

# Load the aircraft 
fdm.load_model('737') 

cutoff_cmd = fdm['propulsion/cutoff_cmd']

# Set engines running
fdm['propulsion/engine[0]/set-running'] = 1
fdm['propulsion/engine[1]/set-running'] = 1

gear_cmd_norm = fdm['gear/gear-cmd-norm']
gear_pos_norm = fdm['gear/gear-pos-norm']

results = []

# Trim
try:
    fdm['simulation/do_simple_trim'] = 2    # Ground trim
    fdm['fcs/throttle-cmd-norm'] = 1.0      # Advance to full throttle

    while fdm.run() and fdm.get_sim_time() <= 120:

        gear_cmd_norm = fdm['gear/gear-cmd-norm']
        gear_pos_norm = fdm['gear/gear-pos-norm']
        elevator_pos_norm = fdm['fcs/elevator-pos-norm']
        
        #if fdm['velocities/vc-kts'] > 100:
        #    fdm['propulsion/cutoff_cmd'] = 1  # Cut-off fuel for engine failure (single engine?)
        
        results.append((fdm.get_sim_time(), 
                        fdm['velocities/vc-kts'], 
                        fdm['propulsion/engine[0]/thrust-lbs'], 
                        fdm['propulsion/engine[0]/n1'], 
                        fdm['aero/alpha-deg'], 
                        fdm['attitude/theta-deg'],
                        fdm['forces/fbz-gear-lbs'],
                        fdm['moments/m-total-lbsft'],
                        fdm['forces/fbx-total-lbs'],
                        fdm['forces/fbx-aero-lbs']))

except RuntimeError as e:
    # The trim cannot succeed. Just make sure that the raised exception
    # is due to the trim failure otherwise rethrow.
    if e.args[0] != 'Trim Failed':
        raise

time, speed, thrust, n1, alpha, theta, fbz_gear, m_moment, fbx_total, fbx_aero = zip(*results)

plot(time, speed, 'Speed vs Time', 'Time (s)', 'Speed (kt)')
plot(time, thrust, 'Thrust vs Time', 'Time (s)', 'Thrust (lbf)')
plot(time, n1, 'N1 vs Time', 'Time (s)', 'N1')
plot(time, alpha, 'Alpha vs Time', 'Time (s)', 'Alpha (deg)')
plot(time, theta, 'Theta vs Time', 'Time (s)', 'Theta (deg)')
plot(time, fbz_gear, 'Gear Z Force vs Time', 'Time (s)', 'Gear Z Force (lbf)')
plot(time, m_moment, 'M Moment vs Time', 'Time (s)', 'M Moment (lbfft)')
plot(time, fbx_total, 'Total X Force vs Time', 'Time (s)', 'Total X Force (lbf)')
plot(time, fbx_aero, 'Aero X Force vs Time', 'Time (s)', 'Aero X Force (lbf)')

Imagine my surprise when I saw this plot of airspeed versus time.

My initial suspicion was that something was happening with the engine, even though I hadn't yet put any code in to cut an engine etc. So added plots for N1 and engine thrust. With N1 staying constant the whole time, and thrust varying with speed as expected. So still a question of what was causing the slowdown, with a weird ~34s period.

So I plotted alpha and theta versus time to get a bit more info on what is happening to the aircraft.

With these periodic short pitch ups to ~4 deg.

Suspecting then that maybe something weird is going on with the gear ground reactions. Looks normal from 0s to ~22s in terms of the normal force decreasing as lift increases. Then a change as the speed starts decreasig after ~22s and then a massive spike around ~30s mark.

Looking at the total force acting on the aircraft, excluding gravity. Again, looks as expected heading towards ~22s mark, decreasing slowly based on thrust decreasing, aerodynamic drag increasing, and the small amount of friction from the gear decreasing slightly. But then surely it doesn't make sense just after that, the **increase **from ~18,000lbf to ~29,000lbf, but the airspeed drops during this period!

So now looking at the aero X force, since the thrust force appears to be as expected. So again from 0s to ~22s as expected as drag increases with increasing airspeed, but then a sudden massive spike to ~5,600lbf of thrust!?

Seems a lot to go wrong for such a simple operation 😉

Just in case there was something weird with my 737, although it's the version from the repo, I tried the F-15, and although not showing the specific periodic cycles above, you can see there is something weird going on.

I guess my first question is whether others can reproduce this, I've included my complete Python program above. And secondly, does anyone know what is causing this?

[bcoconni]

I get the exact same results than you 😄 and could improve things by:

• Calling fdm.run_ic() to initialize the aircraft according to the specified IC
• Pushing the throttle of engine 2 fdm['fcs/throttle-cmd-norm[1]'] = 1.0 to get a symmetrical thrust

--- Sean737.py  2024-08-18 14:03:11.677719842 +0200
+++ Sean737v2.py        2024-08-18 14:05:55.139118136 +0200
@@ -19,6 +19,7 @@
 # Set engines running
 fdm['propulsion/engine[0]/set-running'] = 1
 fdm['propulsion/engine[1]/set-running'] = 1
+fdm.run_ic()
 
 gear_cmd_norm = fdm['gear/gear-cmd-norm']
 gear_pos_norm = fdm['gear/gear-pos-norm']
@@ -29,6 +30,7 @@
 try:
     fdm['simulation/do_simple_trim'] = 2    # Ground trim
     fdm['fcs/throttle-cmd-norm'] = 1.0      # Advance to full throttle
+    fdm['fcs/throttle-cmd-norm[1]'] = 1.0   # Advance to full throttle
 
     while fdm.run() and fdm.get_sim_time() <= 120:
 

Most likely the aircraft was bouncing because it had silly IC. This is what I get with the changes above:

I can't remember how many times I forgot calling FGFDMExec.run_ic() as well. It's a classic 😆

[seanmcleod]

Ah, thanks @bcoconni for confirming, and figuring out the issue 😉

In terms of FGFDMExec.run_ic() I've typically used it in the past, for example as below when plotting AoA vs speed for trim at some altitude etc.

for speed in range(120, 460, 10):
    fdm['ic/h-sl-ft'] = 250
    fdm['ic/vc-kts'] = speed
    fdm['ic/gamma-deg'] = 0

    fdm.run_ic() # Initialize the aircraft with initial conditions

    fdm['simulation/do_simple_trim'] = 1
    results.append((fdm['velocities/vc-kts'], fdm['aero/alpha-deg'], fdm['aero/coefficient/CDde']))

But assumed for my runway case that defaults for things like lat, lon, altitude of all 0 (mid-Atlantic runway) would be fine especially in conjunction with performing a ground trim. My code doesn't explicitly set any initial conditions, but your call to FGFDMExec.run_ic() does the trick, so there is obviously some initialisation of the initial conditions that only happens with the call.

In terms of fcs/throttle-cmd-norm I assumed on a multi-engine aircraft that this would apply to both engines if an index wasn't supplied. I did half wonder if this was really the case since the total X force figures did look a bit low, but there were so many other weird results I didn't look too closely. So the aircraft was probably also going around in circles given the asymmetric thrust!

[bcoconni]

But assumed for my runway case that defaults for things like lat, lon, altitude of all 0 (mid-Atlantic runway) would be fine especially in conjunction with performing a ground trim. My code doesn't explicitly set any initial conditions, but your call to FGFDMExec.run_ic() does the trick, so there is obviously some initialisation of the initial conditions that only happens with the call.

After some further investigations, I remembered that FGFDMExec::RunIC() calls FGPropagate::SetInitialState and FGPropagate::InitializeDerivatives which are doing non trivial tasks. So we cannot assume that JSBSim is initialized to meaningful default values prior to calling FGFDMExec::RunIC().

In terms of fcs/throttle-cmd-norm I assumed on a multi-engine aircraft that this would apply to both engines if an index wasn't supplied.

The class SGPropertyNode assumes that fcs/throttle-cmd-norm and fcs/throttle-cmd-norm[0] are the same property: the lack of index is interpreted as the index [0].