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 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 command 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 😉

[agodemar]

@seanmcleod70 Great example here.

Some comments should be corrected, as they refer to rudder while the code is on roll and ailerons.

Here is the Python code I wrote for this.

# ...

# 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

# ...

# Plot the aileron and rudder commands on the secondary y-axis.
# ...

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

# ...

[seanmcleod70]

@agodemar thanks, I had copy pasted from my original rudder kick example which involved both the rudder and ailerons. I've updated the code to snippet comments above.

[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.

[seanmcleod70]

Assuming theta is constant so thetadot is 0

Remember, thetadot $\dot \theta$ isn't the body pitch rate $q$ in the formula for calculating $\dot \phi$, you're mixing up Euler rate angles and body angular rates.

No because Clp is bigger than Clbeta in this aero (and it seems most, because the turn is coordinated so beta is small)

You started off by saying "neutral rudder", for "THIS situation", so no turn coordination. Also don't confuse the coefficient/derivative value with the value of the moment. For example here is a snippet from the 737.

            <function name="aero/coefficient/Clb">
                <description>Roll_moment_due_to_beta</description>
                <product>
                    <property>aero/qbar-psf</property>
                    <property>metrics/Sw-sqft</property>
                    <property>metrics/bw-ft</property>
                    <property>aero/beta-rad</property>
                    <value>-0.09</value>
                </product>
            </function>

            <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 sure, the magnitude of $C_{l_{p}}$ is larger than $C_{l_{\beta}}$, but look at values for the actual moments $L_p$ and $L_{\beta}$ from my graph above for the 737 moments. I updated the graph in the comment above, used moment names rather than coefficient names etc.

Take a look at the values of the moments at time = 6s. Also no turn coordination, but $\beta$ remains pretty small. With $L_{\beta}$ being ~10% larger in magnitude compared to $L_p$.

Moment	Value
$L_{\beta}$	5,621
$L_p$	-5,011
$L_r$	-5,360

You only included the following columns in your MD-82.csv, "Time, P, Q, R, Phi, Theta, Psi", i.e. no moments to plot and compare.

confirmed to me that Clp is the reason it drifts

Remember $L_p$ is a roll dampening moment, so in the following equation to calculate $\dot \phi$, $L_p$ dampens the body roll rate $p$. With a smaller value for $C_{l_{p}}$ then $p$ in the following equation would be larger, and vice-versa if the coefficient was larger.

The main cause for a negative $\dot \phi$ is the term after $p$, i.e. the combination of $q$, $\phi$, $r$ and $\theta$.

$$ \dot \phi = p + (q \sin \phi + r \cos \phi) \tan \theta$$

Remember, with $p = 0$ that would mean $L_p = 0$, and assuming the same values for $q$, $\phi$, $r$ and $\theta$ you would still have a negative $\dot \phi$, it would in fact be more negative, i.e. it would drift even more.

[Octal450]

I did mix up thetadot and q in my reply, but that was just due to me being tired.

The same applies.

If in a bank with thetadot at 0, q is positive. So it plays no inversion factor. q sin phi also makes negative so phidot is negative when p = 0 and thus it drifts further left when phi = negative.

I didn't understand what you meant about moment vs coefficient. The coefficients build up to the moment no?

Clp IS bigger than Clbeta which is why it drifts left. If I make Clbeta bigger it stops or drifts right. Which is why I said the issue is not always very apparent.

That doesn't mean I'm confusing it, its just not the same in every situation I know. But I'm talking about the situations where the drift occurs.

I'll give you a new output with Cl coefficients. What other properties do you want?

The 737 in JSBSim DONT use it! It has that horrible nonsense "yaw damper" from aeromatic. Except is not a yaw damper!!! It should have a washout on r, otherwise it reacts to steady state r in turn and causes the plane to Slip and thus Clbeta is big!!!! This is not correct and I don't why the aeromatic has this nonsense. It needs washout of r with coefficient roughly 0.5 or so for c1 and then it will not slip in turns and thus Clbeta is smaller. Furthermore I have turn coordinator but no integrator on it so beta is even smaller. This is why I said Clbeta is smaller than Clp.

Please do not use wrong for this, can you direct use the MD-82 flight model on my repo?

J

[Octal450]

The aeromatic and every default model in jsb should have washout added. Otherwise the rudder fights the turns and pushes the plane into a slip uncoordinated turn every time. I have no idea why that was added!

J

[seanmcleod70]

Please do not use wrong for this, can you direct use the MD-82 flight model on my repo?

I'm more than happy to make use of your MD-82 FDM if I can run it via JSBSim stand-alone, i.e. it doesn't require FlightGear Nasal support etc. Then in my example Python scripts above I can simply replace the 737 model with your MD-82 model but still trim it, generate scripted inputs and log and plot parameters etc. and then we're both looking at and using the same flight model, particularly when we want to discuss particular values we see etc.

AIRCRAFT_NAME="737"

I'll give you a new output with Cl coefficients. What other properties do you want?

If I can run your MD-82 FDM in JSBSim stand-alone then you don't need to log any new output, my script already does all the required logging. If it's not feasible to run your FDM in JSBSim stand-alone then please log all the roll moments you have in your FDM plus control positions, in addition to what you included in your previous log.

I didn't understand what you meant about moment vs coefficient. The coefficients build up to the moment no?

JSBSim needs to calculate the relevant forces and moments for each axis in the FDM. Typically FDM authors do that by including functions in the relevant aero axis section, e.g. like the two examples I gave above.

Now the first issue is that the function name implies that it is calculating and returning a coefficient, when in actual fact what the function is calculating and returning is a moment, which happens to make use of a particular coefficient, $C_{l_{\beta}}$ in this first example. So the function name should rather be aero/moment/Lbeta.

            <function name="aero/coefficient/Clb">
                <description>Roll_moment_due_to_beta</description>
                <product>
                    <property>aero/qbar-psf</property>
                    <property>metrics/Sw-sqft</property>
                    <property>metrics/bw-ft</property>
                    <property>aero/beta-rad</property>
                    <value>-0.09</value>
                </product>
            </function>

            <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 in this 737 example we have two fixed coefficients, in general there could be a lookup table etc. based on say $\alpha$ etc. for a coefficient lookup for the relevant moment. In this case the coefficients are $C_{l_{\beta}} = -0.09$ and $C_{l_{p}} = -0.4$.

So as I mentioned, in this 737 example the magnitude of the coefficient $C_{l_{p}}$ is larger than $C_{l_{\beta}}$. But you can't assume that the moment $L_p$ is always going to be larger than the $L_{\beta}$ moment.

It all depends on the other parameters in the moment calculation for each one. Once you exclude the common parameters from both then $L_{\beta}$ depends on $\beta$ and $L_p$ depends on the product of $p$ and aero/bi2vel which is $\frac{b}{2u}$.

So as I pointed out for the 737 example above you had $L_{\beta} = 5,621$ and $L_p = -5,011$, i.e. the magnitude of $L_{\beta}$ was greater than the magnitude of $L_p$.

P.S
Talking of aircraft rolling, you might find the following interesting:

https://seanmcleod70.github.io/FlightDynamicsCalcs/A4SkyhawkRollPerformance.html

https://seanmcleod70.github.io/FlightDynamicsCalcs/FighterRollRates.html

[Octal450]

Hi,
Well it doesn't rely on nasal, but the FCS section relies on some aircraft systems and those may rely on various props initialized in the -set file. I can strip it all out I'm sure or just replace the FCS with a simple one.

OK, I understand your explanation. I also use the coefficient naming scheme. Fair enough.

Regardless, what I am saying is, in the instance of this issue, Clbeta is always less than Clp. If it wasn't, then this issue would not occur. This can be verified by increasing Clbeta coefficient number (0.04 in your 737).

Assuming no other coeffs are at play here (Clr, Clda, Cldr, Cldsp, etc), and the aircraft is not pushed into a skid (beta angle higher) by rudder, or wind, then Clbeta will be trying to roll you level (because the turn coordinator is never perfect at tracking sideslip, so there is always a bit of slipping going on - but NOT skidding unless some other influence like rudder or wind does so) (if there is no turn coordinator, the plane nautrally slips, NOT skids for any yaw-stable airplane), then the roll moment here consists of Clbeta and Clp. Thus in order for this to occur, the magnitude of Clp must be larger than Clbeta and Clp is at fault.

Yes, I DO understand that not every situation it is like this. What I am saying is that in order for what I am describing to occur, |Clp| > |Clbeta|, because the plane is not skidding, there is 0 aileron or rudder input, wind is removed (calm weather), and theta and r are such that the p calculation return a positive (in a left bank) value when phidot is 0.

I am providing you the set of conditions in which the behavior that this issue was about occurs - I of course understand that it is not UNIVERSALLY like this, which is why the issue does not occur universally, it occurs at specific instances which cause p to be positive when phidot is 0, which results in Clp negative, which results in a further left drift, because Clbeta positive is not big enough to overcome the Clp. This is also why offsetting p by phidot or downright replacing the property removes the issue, which isolates Clp and thus p-rad_sec as the cause of this drift.

Again restating - this may be correct behavior. I'm not 100% convinced that this is actually a problem and perhaps it's completely correct to do so, I just thought it was odd.

Btw about the yaw damper thing, I'm going to open a separate issue and I can probably provide patches to add the washout and modify aeromatic in order to it actually be a yaw damper, and not fighting and slipping the plane on every turn - I really don't know how such an oversight was missed for 10+ years but I see aeromatic FDE's all the way back to the early 2010s with this same mistake and it results in very bad behavior. But with such an issue, of course you will NOT see it, because beta is larger than it should be by that bad yaw damper thing (or even with no YD/TC, the issue is less apparent. It is only really apparent in reasonably coordinated turns because beta is small).

I will provide an updated CSV file for you as soon as I am home with many more relevant properties and all moments logged.

Kind Regards,
Josh

[Octal450]

PS: Your write-ups are excellent.

[seanmcleod70]

conditions in which the behavior that this issue was about occurs
But with such an issue, of course you will NOT see it,

Maybe we need to confirm the behavior/it.

My understanding from your original post is that the behavior you're describing is:

• Pitch up after take-off
• Roll left with aileron input
• Zero out aileron input
• Note that body roll rate $p$ is positive while Euler bank angle rate $\dot \phi$ is negative

I reproduced the behavior as per what I've just listed as my understanding of the behavior you were describing by using the 737 model from the JSBSim repo on the 17th Sep above. Here is a link to the graph from that comment to show the 737 matching my understanding of your behavior.

From ~4.5s we end up in a fairly steady state situtation with $p$ being positive and $\dot \phi$ being negative.

[Octal450]

Sean, here is the issue:

Clp here will be negative.

It's a slow but noticeable roll if you pay attention. The AP can also be seen when used slightly with right aileron.

This causes a slow continued left roll in that state:

Pitch up after take-off
Roll left with aileron input
Zero out aileron input
Note that body roll rate p is positive while Euler bank angle rate ϕdot is negative

Kind Regards,
Josh

[seanmcleod70]

Yep, the red box you've drawn over the end portion of my graph pretty much matches what I mentioned:

From ~4.5s we end up in a fairly steady state situtation with $p$ being positive and $\dot \phi$ being negative.

So we're in agreement.

Initially in terms of your post you were surprised that the sign of $p$ and $\dot \phi$ could be opposite. You suspected that a negative $L_p$ roll damping moment was the cause in terms of a positive $p$ resulting in a negative $\dot \phi$.

The graph above is based on the 737 as-is in the repo with a $C_{l_{\beta}} = -0.09$. Watch what happens when I vary $C_{l_{\beta}}$ in terms of the resulting $p$, $\dot \phi$ and how $\dot \phi$ becomes more negative as the magnitude of $L_p$ decreases.

$C_{l_{\beta}}$	$p$	$\dot \phi$	$L_{\beta}$	$L_p$
-0.09	0.44	-0.30	5,650	-4,910
-0.045	0.19	-0.61	4,060	-2,100
-0.0225	0.04	-0.79	2,490	-330
-0.019	0.01	-0.83	2,180	-63

$C_{l_{\beta}} = -0.045$

$C_{l_{\beta}} = -0.019$

[Octal450]

It seems you were right. Well then I have no idea why it rolls left. I guess my knowledge/logic was wrong. Thoughts?

MD-82.csv

J

[agodemar]

@Octal450 @seanmcleod70 Have you guys considered the effect of yaw rate $r$ on roll moment $l$?
If you have a nonzero $r$, the aerodynamics says that the airframe develops an amount of roll moment
$\Delta l (r) = \bar{q}Sb\, C_{lr}\, r$
where the coefficient $C_{lr}$ is usually positive (in general, it depends upon wing aspect ratio, twist, taper ratio, and sweep).
This means that if you're yawing right, the aircraft tends to roll right: the left wing's airfoils experience a greater dynamic pressure than their right wing's counterparts.
From the above plots, I see that the $r(t)$ is mostly negative, hence the $\Delta l(r(t))$ develops into a continuously negative roll moment. That's a roll-to-left tendency that can be significant.

[seanmcleod70]

@agodemar yep, in my plots of the moments for the 737 example above you can see the roll moment due to yaw rate $r$, the green line.

[agodemar]

@seanmcleod70 I see. Yet another evidence of roll-yaw coupling. Well done.

[seanmcleod70]

Well then I have no idea why it rolls left. I guess my knowledge/logic was wrong. Thoughts?

My guess is that you're struggling to visualize things between the two frames, the body frame and the Euler frame.

I tried to give some extreme examples in this comment above to demonstrate how different the body rates and Euler rates can be - #1334 (comment)

Here is another simple example.

Let's assume:

• aircraft is pitched up with $\theta = 45 \mathrm{deg}$, which means $\tan \theta = 1$.
• aircraft has a current Euler bank angle of $\phi = -45 \mathrm{deg}$, which means $\sin \phi = -0.707$
• that the net moments are such that $L = 0$ and the body rates are $p = 0$ and $r = 0$ and that $q = 10 \mathrm{deg/s}$

So given all that, in particular $p = 0$ and the roll moment $L = 0$, would you be surprised that the aircraft's Euler bank rate $\dot \phi$ would be non-zero and negative at $-7.07 \mathrm{deg/s}$, i.e. that the aircraft would increase it's Euler bank angle to the left?

Whereas you seemed to be fixated on trying to tie a single moment, the roll damping moment, in the body frame, to explain why $\dot \phi$ is negative while the body roll rate $p$ is positive.

Rather look at it in terms of the current aircraft state, plus the net moments in all 3 axes of the aircraft in the body frame resulting in a combination of $\theta$, $\phi$ and the body rates $p$, $q$, $r$ and given these 6 states you can see from the formula above how you can end up with a positive body roll rate $p$ and a negative Euler bank rate $\dot \phi$.

If you want to dig a bit deeper for your particular example for a given $\theta$ and $\phi$, then you need to understand why $q$ is positive, and why $r$ is negative, and they happen to be large enough that the sum of the terms that they're part of is greater than $p$.

[seanmcleod70]

@Octal450, @agodemar I noticed a slight difference in what I was expecting the net rolling moment to be versus what I was calculating by adding up all the individual rolling moments, so I took a closer look to figure out the difference and wrote it up here - Moment Accounting.