Calculation of angular accelerations

[matheusdarocha]

Good morning,

I am trying to calculate the angular acceleration using JSBSim output and comparing the answer with the output one.
My inputs:
'velocities/p-rad_sec': -0.008649968842503762
'velocities/q-rad_sec': 0.19870770462160356
'velocities/r-rad_sec': 0.0723236898016243
'moments/l-aero-lbsft': 1895.5918976271505
'moments/m-aero-lbsft': 26534.681031929387
'moments/n-aero-lbsft': -1614.7447248461817

At the JSBSim source code, the formula to calculate the accelerations are:

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

J is the inertia matrix
Jinv is the inverse inertia matrix
vMoments is the moment vector in the body frame
in.vPQRi is the total inertial angular velocity of the vehicle

in.Ti2b is the ECI to body frame transform
in.vOmegaPlanet is the Earth angular rate - expressed in the inertial frame

How can I get in.Ti2b and in.vOmegaPlanet?
and is there a way to remove the angular rate from JSBSim model?

Could you please help me with that?

[matheusdarocha]

[matheusdarocha]

I found that in.Ti2b can be calculated using

q1 = +(math.sin(phi/2.) * math.cos(theta/2.) * math.cos(psi/2.) - math.cos(phi/2.) * math.sin(theta/2.) * math.sin(psi/2.))
q2 = +(math.cos(phi/2.) * math.sin(theta/2.) * math.cos(psi/2.) + math.sin(phi/2.) * math.cos(theta/2.) * math.sin(psi/2.))
q3 = +(math.cos(phi/2.) * math.cos(theta/2.) * math.sin(psi/2.) - math.sin(phi/2.) * math.sin(theta/2.) * math.cos(psi/2.))
q0q0 = q0q0
q1q1 = q1q1
q2q2 = q2q2
q3q3 = q3q3
q0q1 = q0q1
q0q2 = q0q2
q0q3 = q0q3
q1q2 = q1q2
q1q3 = q1q3
q2q3 = q2q3

mT = [[0,0,0],[0,0,0],[0,0,0]]
mT[0][0] = q0q0 + q1q1 - q2q2 - q3q3
mT[0][1] = 2.0*(q1q2 + q0q3)
mT[0][2] = 2.0*(q1q3 - q0q2)
mT[1][0] = 2.0*(q1q2 - q0q3)
mT[1][1] = q0q0 - q1q1 + q2q2 - q3q3
mT[1][2] = 2.0*(q2q3 + q0q1)
mT[2][0] = 2.0*(q1q3 + q0q2)
mT[2][1] = 2.0*(q2q3 - q0q1)
mT[2][2] = q0q0 - q1q1 - q2q2 + q3q3
Ti2b = np.matrix(mT)```

and vOmegaPlanet is [0,0,0.0000729211]

But I can't get the accelations equal to

pp = fdm['accelerations/pdot-rad_sec2'] qp = fdm['accelerations/qdot-rad_sec2'] rp = fdm['accelerations/rdot-rad_sec2']

[seanmcleod]

Are you assuming that in.Moment = moments/l-aero-lbsft + moments/m-aero-lbsft + moments/n-aero-lbsft? What about other moments, for example thrust moments etc.?

[matheusdarocha]

I changed it from aero to total, but they have the same value

# Moments
l = fdm['moments/l-total-lbsft']
m = fdm['moments/m-total-lbsft']
n = fdm['moments/n-total-lbsft']`

[seanmcleod]

If the aero and total moments are identical then either you're using a glider or it so happens your engine is vertically aligned with the cg and for a single engine also horizontally aligned, and for multi-engine they're producing the same thrust.

[seanmcleod]

Also note this distinction.

    /// Angular velocities of the body with respect to the ECI frame (expressed in the body frame).
    FGColumnVector3 vPQRi;
    /// Angular velocities of the body with respect to the local frame (expressed in the body frame).
    FGColumnVector3 vPQR;

  PropertyManager->Tie("velocities/p-rad_sec", this, eP, (PMF)&FGPropagate::GetPQR);
  PropertyManager->Tie("velocities/q-rad_sec", this, eQ, (PMF)&FGPropagate::GetPQR);
  PropertyManager->Tie("velocities/r-rad_sec", this, eR, (PMF)&FGPropagate::GetPQR);

  PropertyManager->Tie("velocities/pi-rad_sec", this, eP, (PMF)&FGPropagate::GetPQRi);
  PropertyManager->Tie("velocities/qi-rad_sec", this, eQ, (PMF)&FGPropagate::GetPQRi);
  PropertyManager->Tie("velocities/ri-rad_sec", this, eR, (PMF)&FGPropagate::GetPQRi);

[matheusdarocha]

The code I am using is showed below if it helps. I am using the ECI to calculate vPQRi

# Calculate the ECI to body frame rotation matrix
def f_Ti2b(phi, theta, psi):
    # phi is the roll attitude angle
    # theta is the pitch attitude angle
    # psi is the yaw attitude angle
    
    # According to JSBSim Source code, file: FGQuaternion.cpp Line 106
    # void FGQuaternion::InitializeFromEulerAngles(double phi, double tht, double psi)
    q0 = +(math.cos(phi/2.) * math.cos(theta/2.) * math.cos(psi/2.) + math.sin(phi/2.) * math.sin(theta/2.) * math.sin(psi/2.))
    q1 = +(math.sin(phi/2.) * math.cos(theta/2.) * math.cos(psi/2.) - math.cos(phi/2.) * math.sin(theta/2.) * math.sin(psi/2.))
    q2 = +(math.cos(phi/2.) * math.sin(theta/2.) * math.cos(psi/2.) + math.sin(phi/2.) * math.cos(theta/2.) * math.sin(psi/2.))
    q3 = +(math.cos(phi/2.) * math.cos(theta/2.) * math.sin(psi/2.) - math.sin(phi/2.) * math.sin(theta/2.) * math.cos(psi/2.))
    
    # According to JSBSim Source code, file: FGQuaternion.cpp Line 187
    # void FGQuaternion::ComputeDerivedUnconditional(void) const
    q0q0 = q0*q0
    q1q1 = q1*q1
    q2q2 = q2*q2
    q3q3 = q3*q3
    q0q1 = q0*q1
    q0q2 = q0*q2
    q0q3 = q0*q3
    q1q2 = q1*q2
    q1q3 = q1*q3
    q2q3 = q2*q3

    mT = [[0,0,0],[0,0,0],[0,0,0]]
    mT[0][0] = q0q0 + q1q1 - q2q2 - q3q3
    mT[0][1] = 2.0*(q1q2 + q0q3)        
    mT[0][2] = 2.0*(q1q3 - q0q2)
    mT[1][0] = 2.0*(q1q2 - q0q3)
    mT[1][1] = q0q0 - q1q1 + q2q2 - q3q3
    mT[1][2] = 2.0*(q2q3 + q0q1)
    mT[2][0] = 2.0*(q1q3 + q0q2)
    mT[2][1] = 2.0*(q2q3 - q0q1)
    mT[2][2] = q0q0 - q1q1 - q2q2 + q3q3
    Ti2b = np.matrix(mT)
    return Ti2b

# Inertial Moments
Ixx = fdm['inertia/ixx-slugs_ft2']
Iyy = fdm['inertia/iyy-slugs_ft2']
Izz = fdm['inertia/izz-slugs_ft2']
Ixy = fdm['inertia/ixy-slugs_ft2']
Ixz = fdm['inertia/ixz-slugs_ft2']
Iyz = fdm['inertia/iyz-slugs_ft2']

# Inertial Matrix
J = np.matrix([[Ixx, -Ixy, -Ixz], [-Ixy, Iyy, -Iyz], [-Ixz, -Iyz, Izz]])
Jinv = np.linalg.inv(J)

# Velocities at the body frame
p = fdm['velocities/p-rad_sec']
q = fdm['velocities/q-rad_sec']
r = fdm['velocities/r-rad_sec']

# Velocities at the ECI frame
pi = fdm['velocities/pi-rad_sec']
qi = fdm['velocities/qi-rad_sec']
ri = fdm['velocities/ri-rad_sec']
vPQRi = np.matrix([p, q, r]).transpose()  # rad_sec

# According to eq 1.3-9 page 22 Stevens and Lewis, 
# "Aircraft Control and Simulation" Second edition (2004)
omega = np.matrix([[0, -r, q], [r, 0, -p], [-q, p, 0]]) # rad / sec

# Moments
l = fdm['moments/l-total-lbsft']
m = fdm['moments/m-total-lbsft']
n = fdm['moments/n-total-lbsft']
Moment = np.matrix([l, m, n]).transpose() # sl * ft / sec²

# According to JSBSim Source code, file: FGAcceleration.cpp Line 161
# vPQRidot = in.Jinv * (in.Moment - in.vPQRi * (in.J * in.vPQRi));
vPQRidot = Jinv * (Moment - (omega*(J*vPQRi)))          # slug * ft² / sec² * slug * ft²


# earth rotation rate at ECI frame
# According to JSBSim Source code, file: FGInertial.cpp 
RotationRate    = 0.00007292115 # Line 56
vOmegaPlanet = np.matrix([0, 0, RotationRate]).transpose() # Line 72

# attitude angles
phi = fdm['attitude/phi-rad']
theta = fdm['attitude/theta-rad']
psi = fdm['attitude/psi-rad']

# ECI to body frame transform
Ti2b = f_Ti2b(phi, theta, psi)

# According to JSBSim Source code, file: FGAcceleration.cpp Line 162
# vPQRdot = vPQRidot - in.vPQRi * (in.Ti2b * in.vOmegaPlanet);
vPQRdot = vPQRidot - omega * (Ti2b * vOmegaPlanet);

# JSBSim output accelerations
pp = fdm['accelerations/pdot-rad_sec2']
qp = fdm['accelerations/qdot-rad_sec2']
rp = fdm['accelerations/rdot-rad_sec2']
PQR_JSBSim = np.matrix([pp, qp, rp]).transpose()

print(f"JSBSim output acceleration:\n{PQR_JSBSim}\n\nCalculated Output acceleration:\n{vPQRdot}")

The output found is showed below

JSBSim output acceleration:
[[3.02134033e-05]
 [3.67735802e-06]
 [1.36766898e-06]]

Calculated Output acceleration:
[[-0.00114528]
 [ 0.00061324]
 [-0.00149393]]

[matheusdarocha]

I was doing some tests and discovered that the follow variables changes A LOT the value of the angular accelerations

fdm['simulation/integrator/rate/rotational'
fdm['simulation/integrator/rate/translational'
fdm['simulation/integrator/position/rotational'
fdm['simulation/integrator/position/translational'

But, looking at the source code, they just change the integration method.
In my opinion, the integration method should not change the value of the angular accelerations. Am I wrong?

[seanmcleod]

In my opinion, the integration method should not change the value of the angular accelerations. Am I wrong?

vPQRidot = in.Jinv * (in.Moment - in.vPQRi * (in.J * in.vPQRi));

So the angular accelerations computed use the angular velocities as input. Which means different integration schemes may result in a different angular velocity for the following time step which means the angular accelerations will differ in the subsequent time step.

[bcoconni]

How can I get in.Ti2b and in.vOmegaPlanet?

The matrices Ti2b and Tb2i can be accessed by FGPropagate::GetTi2b() and FGPropagate::GetTb2i() if you're using C++.
The vector vOmegaPlanet is supplied by FGInertial::GetOmegaPlanet().

Neither of them are accessible using Python unfortunately but this can be added quite easily.

[bcoconni]

and is there a way to remove the angular rate from JSBSim model?

If you are talking about nullifying vOmegaRate you can add a <planet> XML element to your FDM:

<planet>
  <rotation_rate> 0.0 </rotation_rate>
</planet>