P-rad_sec behavior

[Octal450]

Hi,
I've noticed this for a long time (many years) but finally wondering about it. At slower speeds mainly, I noticed that when attitude is stable, the roll rate reported by FlightGear and JSBsim disagree. At higher speeds the issue doesn't occur and the two values agree much much closer.

The most extreme case here just after takeoff in a 20 degree LEFT bank, where the signs differ, FlightGear reports the aircraft is very slowly rolling further left (which it is), but JSBsim p-rad_sec (and aero) reports the plane is slowly rolling right. This causes Clp coefficient to be strongly negative which is causing the plane to keep rolling left...

This means the planes are roll attitude unstable just slightly, when they should really be attitude stable or even rolling slowly back to center. Just Clp is wrong here, Clr is correctly negative here in the bank. Clbeta is trying to roll right as it should. Surface moments are tiny as the controls are centered and the yaw damper is off.

What is going on here? Shouldn't p-rad_sec (and aero) be not disagreeing like this? If the plane is rolling slowly left, then p-rad_sec (and aero) should be showing a left roll. It's not a polarity issue, with larger banks, its close enough.

When I replace p-aero-rad_sec with /orientation/roll-rate-degps in the aero Clp, it behaves as I expect (minus missing wind info)

Kind Regards,
Josh

[seanmcleod70]

Hi Josh

What is going on here? Shouldn't p-rad_sec be not disagreeing like this?

The issue is that p-rad_sec and roll-rate-degps aren’t the same thing .

The roll-rate-degps is a Euler angular rate and p-rad_sec is a body angular rate.

Taken from FlightGear GitLab:

FGJSBsim::copy_from_JSBsim()

    _set_Omega_Body( Propagate->GetPQR(FGJSBBase::eP),
                     Propagate->GetPQR(FGJSBBase::eQ),
                     Propagate->GetPQR(FGJSBBase::eR) );

    _set_Euler_Rates( Auxiliary->GetEulerRates(FGJSBBase::ePhi),
                      Auxiliary->GetEulerRates(FGJSBBase::eTht),
                      Auxiliary->GetEulerRates(FGJSBBase::ePsi) );

void FlightProperties::set_Euler_Rates(double x, double y, double z)
{
  _root->setDoubleValue("orientation/roll-rate-degps", x * SG_RADIANS_TO_DEGREES);
  _root->setDoubleValue("orientation/pitch-rate-degps", y * SG_RADIANS_TO_DEGREES);
  _root->setDoubleValue("orientation/yaw-rate-degps", z * SG_RADIANS_TO_DEGREES);
}

If you’re not familiar with the difference take a look at - https://aviation.stackexchange.com/questions/83993/the-relation-between-euler-angle-rate-and-body-axis-rates

Cheers

[Octal450]

Hi Sean,
Yes, but my question is why is it causing this indication?

For Q I understand, it behaves how I expect. Because in banks the lift. So the 0Q is negative Euler rate.

Otherwise I don't know how to solve this aero issue. Cause at lower speeds it banks away from center when at bigger bank angles due to Clp. Aside from just arbitrarily increasing Clbeta or reducing Clr. Which doesn't seem right. It's not just my planes either.

Kind Regards,
Josh

[bcoconni]

Yes, but my question is why is it causing this indication?

As @seanmcleod70 explained: they are not the same thing. If you read the text at the link he has provided, you will find the relation between the parameters:

The property p-rad_sec is $p$ in the matrix equation above and roll-rate-degps is $\dot{\phi}$ and you have the relationship:

$$p=\dot{\phi}-\sin\theta\dot{\psi}$$

Hence the 2 properties are different.

Otherwise I don't know how to solve this aero issue. Cause at lower speeds it banks away from center when at bigger bank angles due to Clp. Aside from just arbitrarily increasing Clbeta or reducing Clr. Which doesn't seem right.

Given that p-rad_sec is 0.00447 and Clp is -8899.97, this means that there is a ratio of −1991045 between the 2 values. So unless you are using a coefficient that big, I'm struggling to connect the dots the way you do.

We definitely need to understand how Clp is computed before jumping to a conclusion.

[seanmcleod70]

Josh can you:

• Provide a link to your FDM
• 5sec graph with (bank angle, p, r, beta, individual Cl moments, Cl total)

[Octal450]

I'll get that to you.

https://github.com/Octal450/MD-80/blob/master/FDE/Config/md-80-aerodynamics.xml#L531

NOTE: this issue is NOT isolated to my planes, and occurs in most JSBsim aircraft in FlightGear that are plausible.

I can get the graphs to you but not this second.

Bests,
Josh

[bcoconni]

I'll get that to you.

Thanks.

The first thing that strikes me is the usage of the property aero/bi2vel:

<function name="aero/coefficient/Clp">
	<description>Roll moment due to roll rate</description>
	<product>
		<property>aero/qbar-psf</property>
		<property>metrics/Sw-sqft</property>
		<property>metrics/bw-ft</property>
		<property>aero/bi2vel</property>
		<property>velocities/p-aero-rad_sec</property>
		<value>-0.51</value>
	</product>
</function>

This property is basically the ratio of a constant (the wing span) to the total velocity. I would expect this value to go through the roof at extremely low speed. Although this is supposed to be compensated by the multiplication with aero/qbar-psf which is proportional to $V^2$. However weird things can happen with finite precision arithmetic.

jsbsim/src/models/FGAerodynamics.cpp

Line 141 in 405c00d

const double twovel=2*in.Vt;

jsbsim/src/models/FGAerodynamics.cpp

Lines 158 to 161 in 405c00d

	if (twovel != 0) {
	bi2vel = in.Wingspan / twovel;
	ci2vel = in.Wingchord / twovel;
	}

In addition to the parameters requested by @seanmcleod70, could you also give the values for all the parameters that are used to compute aero/coefficient/Clp ?

[bcoconni]

And indeed the constant is quite large. The moment is

$$M_{l_p}=\dfrac{1}{2}\rho SV^2\dfrac{\bar{b}^2}{V}pC_{l_p}$$

which gives $Sb^2=14057416\text{ ft}^4$ with the values from the MD80 ( $S=1209\text{ ft}^2$ and $b=107.83\text{ ft}$ )

[seanmcleod70]

This property is basically the ratio of a constant (the wing span) to the total velocity. I would expect this value to go through the roof at extremely low speed.

I doubt that's an issue here since I'm pretty sure Josh is referring to his airliner model MD-80 soon after take-off, so I would imagine around the 150KIAS range or so. As opposed to the issue we had a while ago with a stationary aircraft and wind and things like alpha-dot etc. blowing up.

[seanmcleod70]

So I've gone ahead and used the 737 in the JSBSim repo to model a roll to a particular bank angle at fairly low speed and plotted the results of the various parameters I was asking Josh for from his example.

At higher speeds the issue doesn't occur...

At lower speeds we're more than likely talking about a higher AoA and assuming no auto-coordination during the bank/turn it means some of the alpha is being converted into beta.

The most extreme case here just after takeoff in a 20 degree LEFT bank

So in my example I've trimmed the 737 at 160KIAS (no flaps) which results in an AoA of 11.5 degrees.

With regards to your initial comment about a sign difference between velocities/p-aero-rad_sec and velocities/phidot-rad_sec you'll see that for 0.5s between 4.5s and 5.0s.

So it will be interesting to see how different your MD-80 example may or may not be.

Here is the Python code I wrote for this.

import jsbsim
import matplotlib.pyplot as plt
import math

PATH_TO_JSBSIM_FILES="../JSBSim"

AIRCRAFT_NAME="737"

# Avoid flooding the console with log messages
#jsbsim.FGJSBBase().debug_lvl = 0

# Create a flight dynamics model (FDM) instance.
fdm = jsbsim.FGFDMExec(PATH_TO_JSBSIM_FILES)

# Load the aircraft model
fdm.load_model(AIRCRAFT_NAME)

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

# Set alpha range for trim solutions
fdm['aero/alpha-max-rad'] = math.radians(12)   # Maximum angle of attack in radians.
fdm['aero/alpha-min-rad'] = math.radians(-4.0) # Minimum angle of attack in radians.

# Initial conditions
fdm['ic/h-sl-ft'] = 200   # altitude above sea level (ft)
fdm['ic/vc-kts'] = 160    # calibrated airspeed (kts)
fdm['ic/gamma-deg'] = 0   # flight path angle (deg)

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

# Attempt to trim the aircraft.
try:
    fdm['simulation/do_simple_trim'] = 1
except jsbsim.TrimFailureError:
    print("Trim failed, continuing rudder kick in an untrimmed state.")
    pass  # Ignore trim failure

times = []      # List to record the simulation time at each step.
betas = []      # List to record the beta angle at each step.
alphas = []
bankAngle = []  # List to record the bank angle at each step.

ps = []
phidots = []

Clbs = []
Clps = []
Clrs = []
Cldas = []
ClTotals = []

ailerons = []   # List to record the aileron control surface deflection.

fdm['fcs/aileron-cmd-norm'] = -0.3

while(fdm.get_sim_time() < 8):
    fdm.run()

    # Record the simulation data.
    times.append(fdm.get_sim_time())
    betas.append(fdm['aero/beta-deg'])
    alphas.append(fdm['aero/alpha-deg'])
    bankAngle.append(fdm['attitude/phi-deg'])
    ailerons.append(fdm['fcs/aileron-cmd-norm'])
    ps.append(math.degrees(fdm['velocities/p-aero-rad_sec']))
    phidots.append(math.degrees(fdm['velocities/phidot-rad_sec']))

    Clb = fdm['aero/coefficient/Clb']
    Clp = fdm['aero/coefficient/Clp']
    Clr = fdm['aero/coefficient/Clr']
    Clda = fdm['aero/coefficient/Clda']
    ClTotal = Clb + Clp + Clr + Clda
    
    Clbs.append(Clb)
    Clps.append(Clp)
    Clrs.append(Clr)
    Cldas.append(Clda)
    ClTotals.append(ClTotal)
    
    if fdm['attitude/phi-deg'] < -17:
        fdm['fcs/aileron-cmd-norm'] = 0

plt.rcParams["figure.figsize"] = (12, 8) # Set the figure size.
fig, (ax1, ax3, ax4) = plt.subplots(3)

# Create a secondary y-axis for the control surface positions.
ax2 = ax1.twinx()

# Plot the beta data on the primary y-axis.
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Angles (deg)')
line1 = ax1.plot(times, betas, label='Beta', color='red')
line2 = ax1.plot(times, bankAngle, label='Bank Angle')
#line3 = ax1.plot(times, alphas, label='Alpha')

# Plot the aileron and rudder commands on the secondary y-axis.
ax2.set_ylabel('Inceptor Position')
line4 = ax2.plot(times, ailerons, label='Aileron', color='grey')

# Add a legend to the plot.
ax1.legend(handles=line1+line2+line4, loc='center right')

ax3.plot(times, ps, label='p-aero')
ax3.plot(times, phidots, label='phidot')
ax3.set_ylabel('Rates (deg/s)')
ax3.set_xlabel('Time (s)')
ax3.legend()
ax3.axhline(y=0, color='grey', linestyle='--')

ax4.plot(times, Clbs, label='$Cl_{\\beta}$')
ax4.plot(times, Clps, label='$Cl_p$')
ax4.plot(times, Clrs, label='$Cl_r$')
ax4.plot(times, Cldas, label='$Cl_{\\delta a}$')
ax4.plot(times, ClTotals, label='Cl Total')
ax4.axhline(y=0, color='grey', linestyle='--')
ax4.legend()

# Set the title of the plot.
plt.title('Roll Test')

fig.tight_layout()

# Display the plot.
plt.show()

[bcoconni]

Just for the record, the 737 is using the same computation than @Octal450's MD80 for Clp:

jsbsim/aircraft/737/737.xml

Lines 706 to 716 in ac3ec7c

	<function name="aero/coefficient/Clp">
	<description>Roll_moment_due_to_roll_rate</description>
	<product>
	<property>aero/qbar-psf</property>
	<property>metrics/Sw-sqft</property>
	<property>metrics/bw-ft</property>
	<property>aero/bi2vel</property>
	<property>velocities/p-aero-rad_sec</property>
	<value>-0.4</value>
	</product>
	</function>

So the comparison is fair 😉

[seanmcleod70]

So, at first glance I don't see anything obviously wrong with the 737 response.

The initial aileron input produces a large $Cl_{\delta a}$ moment, which then gets the roll started to the left, and as the roll rate $p$ increases the roll damping $C_{l_{p}}$ then starts increasing to oppose the roll.

At the same time $\beta$ starts increasing. The 737 model doesn't have an adverse/proverse yaw coefficient, $C_{n_{\delta a}}$ so the $\beta$ isn't from that, it's from rolling with a non-zero $\alpha$ and no turn-coordination. I see the MD-80 model also doesn't have $C_{n_{\delta a}}$ defined. Given the increasing $\beta$ there is then an increase in $C_{l_{\beta}}$ which also opposes the roll. There will be a moment via $C_{n_{\beta}}$ trying to reduce $\beta$.

Given the yaw rate increase there is an increase in $C_{l_{r}}$ which contributes to the roll.

Also given the relationship between $p$ and $\dot{\psi}$:

$$p=\dot{\phi}-\sin\theta\dot{\psi}$$

You can then see difference between $p$ and $\dot{\psi}$ in the graph starting round about the 1s mark.

When the aileron is centered as a step response at 3.5s the total $Cl$ moment switches sign so the roll rate $p$ starts slowing and therefore $C_{l_{p}}$ starts diminishing. There is still a small $\beta$ so $C_{l_{\beta}}$ is still generating a moment in the opposite direction.

At 4.5s the roll rate $p$ becomes 0 and therefore $C_{l_{p}}$ also equals 0 at this instant, with the net moment being the sum of $C_{l_{\beta}}$ and $C_{l_{r}}$ and since $C_{l_{\beta}}$ is larger the aircraft now starts rolling back in the opposite direction to the initial roll, but with a fairly small magnitude.

[bcoconni]

In addition, the constant $Sb^2$ is quite similar in amplitude for the 737 and the MD80:

• MD80: $Sb^2=14057416\text{ ft}^4$
• 737: $Sb^2=10501633\text{ ft}^4$

So both aircraft should give a similar roll damping value in the same flying conditions (i.e. same velocity $V$, same air density $\rho$ and same rolling rate $p$).

[seanmcleod70]

Just to be clear, I've been a bit fast and loose in terms of often referencing the moment by it's derivative/coefficient. So when I say $C_{l_{p}}$ increases for example I don't mean the actual derivative/coefficient increases, it's fixed in the 737 model, rather I'm referring to the moment.

The function name is typically badly named in JSBSim FDMs, the function doesn't calculate and return $C_{l_{p}}$ as per the coefficient build-up technique, rather it's using that derivative to calculate a moment, as per the description.

            <function name="aero/coefficient/Clp">
                <description>Roll_moment_due_to_roll_rate</description>
                <product>
                    <property>aero/qbar-psf</property>
                    <property>metrics/Sw-sqft</property>
                    <property>metrics/bw-ft</property>
                    <property>aero/bi2vel</property>
                    <property>velocities/p-aero-rad_sec</property>
                    <value>-0.4</value>
                </product>
            </function>

And lastly derivatives involving the body rates, $p$, $q$, $r$ are generally expressed as non-dimensional rates.

https://www.aircraftflightmechanics.com/StaticStability/Nondimensionalrates.html

So the JSBSim function above really involves $C_{l_{\tilde p}}$ as opposed to $C_{l_{p}}$. Might have to squint to see the small tilde over the $p$.

[bcoconni]

I've run a check based on @seanmcleod70 Python script to verify that $p$ and $\dot{\phi}-\dot{\psi}\sin\theta$ are matching and they are:

Just added the collection of $\theta$ and $ \dot{\psi}$ and the computation of the formula to the script:

    psidots.append(fdm['velocities/psidot-rad_sec'])
    thetas.append(fdm['attitude/theta-rad'])

ps_theory = np.array(phidots)-np.sin(np.array(thetas))*np.array(psidots)

[seanmcleod70]

The first thing I did after Josh shared a link to his FDM was to confirm the set of rolling moments to see if there were any that the 737 model doesn't have compared to the MD80, and then to compare the formula for each to see that the formulas for calculating the individual moments were the same, and then lastly to compare the sign and magnitude of the derivatives/coefficients.

Where necessary I've converted the MD80 derivatives from per deg to per radian.

Derivative/Coefficient	737	MD80	MD80 Comments
$C_{l_{\beta}}$	-0.09	-0.09 to -0.11	For $\alpha$ from -4.4 deg to +6.1 deg
$C_{l_{p}}$	-0.4	-0.51
$C_{l_{r}}$	0.09	0.1

So, I'd be surprised if the MD80's response isn't very similar to the 737's response.

[Octal450]

Forgive me for not being as well versed in aero as you guys are.

My end goal here is that when stabilized at a bank angle, Clp does not "try" to roll the airplane away from neutral by the cause of P - in real I can't see how there would be a moment from that, I'd think the moment would trend towards level wings, not away.

I should note that with no turn coordination enabled, the issue is less apparent because the beta angle causes it to try to roll level. But with turn coordination turned on, it is apparent. IMPORTANT NOTE: The turn coordinator does NOT use an integrator and is not very tight, even in real, so this issue is NOT caused by too much rudder from turn coordination causing a roll moment. I have verified this already - it was my first guess as what could be happening.

I should note that I temporarily put

	<function name="velocities/p-fixed-aero-rad_sec">
		<description>Corrected p-aero-rad_sec</description>
		<sum>
			<toradians><property>/orientation/roll-rate-degps</property></toradians>
			<difference>
				<property>velocities/p-aero-rad_sec</property>
				<property>velocities/p-rad_sec</property>
			</difference>
		</sum>
	</function>

I guess roll-rate-degps is just phidot from JSBsim?

And fed that into Clp and it resulted in expected behavior - using the aero and non aero diff to add wind effects.

Kind Regards,
Josh

[Octal450]

	<wingarea unit="FT2">1209</wingarea>
	<wingspan unit="FT">107.83</wingspan>
	<chord unit="FT">13.39</chord>

Variables used in aero

[seanmcleod70]

I guess roll-rate-degps is just phidot from JSBsim?

No need to guess, I stated that in my very first response 😉

FGJSBsim::copy_from_JSBsim()
{
    _set_Euler_Rates( Auxiliary->GetEulerRates(FGJSBBase::ePhi),
 ...

void FlightProperties::set_Euler_Rates(double x, double y, double z)
{
  _root->setDoubleValue("orientation/roll-rate-degps", x * SG_RADIANS_TO_DEGREES);

Generate a plot like I did for the 737 and let's take a look at the data.

In the meantime can you spot anything you're not expecting in the 737 plot above?

[seanmcleod70]

"A picture tells a thousand words...", but while we wait for the equivalent graph data for the MD80 here is a potential match of the original figures mentioned compared to the 737's graph data.

Parameter	MD80	737 at 4.8s
$\dot \phi$	-0.14	-0.15
$p$	0.26	0.40
$C_{l_{\beta}}$	17,506	15,200
$C_{l_{p}}$	-8,899	-4,500
$C_{l_{r}}$	-8,146	-7,600
Net $C_l$	461	3,100

So they match in terms of their sign in particular, and the magnitudes are fairly similar.

[Octal450]

Hi Sean,
In your plot, the following situation:

Shows a situation where Clp will have the wrong polarity from what I'd expect.

I will do a log for you.

Kind Regards,
Josh

[seanmcleod70]

But $C_{l_{p}}$ doesn't have the wrong polarity here in the 737 graph. When $p$, the blue line is negative $C_{l_{p}}$ is positive and vice-versa, i.e. always opposite since it's a damping term.

[Octal450]

Angle bank is phi and phidot is amber line. So phidot is having opposite sign and different magnitude as p. So Clp is not helping there because of the sign.

Anyways, I've found out, apparently, airliners sometimes at higher pitch angle do this... so maybe I just leave it.

[Octal450]

Hi Sean,
Please see the following:

I had banked to -20 phi on the left and then let it sit there, it drifted slowly over this graph to -25 phi where I then banked it back to -20 phi.

Kind Regards,
Josh