Moment Accounting

[seanmcleod70]

In a recent discussion - P-rad_sec behavior I logged all the roll moments defined in the aircraft's FDM to try and understand the contributions from each and to look at the net rolling moment.

Given that at first glance it looked like the body roll rate $p$ at the end of the trace was fairly close to being constant I assumed that the net total roll moment would also be fairly close to 0.

However on inspection the net roll moment at t = 6s was $L_{\mathrm{total}} = -5,923$, which was larger in magnitude compared to what I was expecting. The 737's $I_{xx} = 562,000$, so that would equate to a body roll rate acceleration of $-0.6 \frac{deg}{s^2}$, which also seemed excessive based on the slope of $p$.

Taking a look at accelerations/pdot-rad_sec2 at t = 6s, showed a value of $-0.16 \frac{deg}{s^2}$, so roughly a 4x difference.

So it did appear that the net rolling moment that I was calculating by adding up all the individual roll moments was too large.

Taking a look at moments/l-aero-lbsft showed a value of -1,875 compared to -5,923 that I had calculated, also moments/l-total-lbsft also showed a value of -1,875, i.e. there weren't other non-aerodynamic moments, e.g. from engines contributing to the final roll moment.

The reason for the difference in net rolling moment is due to the MRC (Moment Reference Center)/AeroRP not being coincident with the current cg.

Which means that the final aerodynamic rolling moment includes an additional moment contribution based on the aerodynamic forces and the moment arm between the current cg and the MRC, as per:

jsbsim/src/models/FGAerodynamics.cpp

Line 288 in 274f1fb

vMoments = vMomentsMRCBodyXYZ + vDXYZcg*vForces; // M = r X F

To double-check I changed the MRC to match the current cg (based on fuel etc.) and now my calculated net rolling moment and moments/l-total-lbsft are identical at 405.

In particular with the 737 model is there is an approximately 4ft difference in the Z-axis between the MRC and the cg, plus this side-force which will result in an additional rolling moment.

jsbsim/aircraft/737/737.xml

Lines 633 to 643 in 274f1fb

	<axis name="SIDE">
	<function name="aero/coefficient/CYb">
	<description>Side_force_due_to_beta</description>
	<product>
	<property>aero/qbar-psf</property>
	<property>metrics/Sw-sqft</property>
	<property>aero/beta-rad</property>
	<value>-1</value>
	</product>
	</function>
	</axis>

So something to watch out for if your moments you're logging and comparing don't add up like you think they should 😉

[seanmcleod70]

In terms of accounting for roll acceleration $\dot p$, dividing the net moment by $I_{xx}$ like I did above is only a first order approximation.

Note the additional term, the cross-product of the angular velocity with the aircraft's angular momentum in the rotational dynamic equation:

$$\mathbf{M}_B = [I_B] \frac{{}^B d\underline{\omega}_B}{dt} + \underline{\omega}_B \times [I_B] \underline{\omega}_B$$

$$\underline{\omega}_B=(p, q, r)$$

$$\frac{d\underline{\omega}_B}{dt} = [I_B]^{-1} ( \mathbf{M}_B -[\underline{\tilde\omega}_B] [I_B] \underline{\omega}_B )$$

JSBSim implements this here:

jsbsim/src/models/FGAccelerations.cpp

Lines 162 to 163 in 274f1fb

	vPQRidot = in.Jinv * (in.Moment - in.vPQRi * (in.J * in.vPQRi));
	vPQRdot = vPQRidot + in.vPQRi * (in.Ti2b * in.vOmegaPlanet);

Assuming a left-right symmetric aircraft (most are) then a number of entries in the inertia tensor matrix disappear, resulting in the following system of equations:

$$\mathcal{L} = \dot p I_{xx} + qr (I_{zz} –I_{yy}) –(\dot r + pq) I_{xz}$$

$$\mathcal{M} = \dot q I_{yy} –pr (I_{zz} –I_{xx}) + (p^2 + r^2) I_{xz}$$

$$\mathcal{N} = \dot r I_{zz} + pq (I_{yy}-I_{xx}) + (qr-\dot p) I_{xz}$$

Solving for the angular accelerations:

$$\dot p = \frac{I_{zz}\mathcal{L} + I_{xz}\mathcal{N} –{ I_{xz}(I_{yy}-I_{xx}-I_{zz})p + [I_{xz}^2 + I_{zz}(I_{zz}-I_{yy})]r}q}{I_{xx}I_{zz}-I_{xz}^2}$$

$$\dot q = \frac{\mathcal{M} –(I_{xx}-I_{zz})pr –I_{xz})p^2-r^2)}{I_{yy}}$$

$$\dot r = \frac{I_{xz}\mathcal{L}+I_{xx}\mathcal{N}+{I_{xz}(I_{yy}-I_{xx}-I_{zz})r+[I_{xz}^2+I_{xx}(I_{xx}-I_{yy})]p}q}{I_{xx}I_{zz}-I_{xz}^2}$$

The $I_{xz}$ component is typically not zero, but is usually smaller than the other components, so if it is neglected then you get the simplified form of Euler's equations:

$$\dot p = \frac{\mathcal{L} + (I_{zz} –I_{yy})qr}{I_{xx}}$$

$$\dot q = \frac{\mathcal{M} + (I_{zz} –I_{xx})pr}{I_{yy}}$$

$$\dot r = \frac{\mathcal{N} + (I_{xx} –I_{yy})qp}{I_{zz}}$$

I used Flight Stability and Automatic Control by Nelson which has a detailed derivation of the rigid body equations of motion.

The following website has a summarized version, but usefully includes the relevant equations in Latex form, https://academicflight.com/articles/equations-of-motion/