P6DOF Theory

Introduction

This document exists to explain the physics of the P6DOF Mover (see the WSF_P6DOF_MOVER input type and the WsfP6DOF_Mover script class). The kinematic and dynamic equations are defined, and where applicable, parallels are drawn between the mathematical symbols in this document and existing P6DOF commands.

The WSF_P6DOF_MOVER provides a 6 DOF (degrees-of-freedom) physics-based mover that supports a wide range of aircraft and weapons modeling. The 6DOFs include translational position (\vec{x}) as well as attitude (\psi, \theta, \phi). The attitude of the vehicle is determined by momentum-based calculations, providing realistic vehicle pointing as well as proper angle of attack (\alpha) and angle of sideslip (\beta) modeling. Rotational kinematics are influenced by aerodynamic and propulsion forces and moments along with rotational inertia. Although P6DOF movers include full 6DOF modeling, they utilize some angular/attitude simplifications to reduce the required data and to simplify control – hence the “pseudo” of the Pseudo-6DOF name. In particular, they only use the diagonal of the inertia tensor to reduce inertial cross-coupling (I_{xx}, I_{yy}, I_{zz}). However, P6DOF movers support a full 6 DOFs and include a detailed, stability derivative build-up approach to aerodynamics as described in this document. This allows detailed aerodynamics modeling, based on Mach, to provide accurate modeling of the transonic, supersonic, and hypersonic regimes.

P6DOF movers include type-specific modeling of turbojets, turbofans, ramjets/scramjets, liquid-propellant rockets, and solid-propellant rockets, providing a wide range of propulsion types. This means each engine model includes the parameters that have the largest impact on thrust performance, such as Mach and altitude.

List of Symbols

Symbol

Definition

P6DOF Command

P6DOF Script Method

\alpha

angle of attack

GetAlpha

\beta

angle of side slip

GetBeta

M

Mach number

GetMach

\psi

yaw angle

GetHeading

\theta

pitch angle

GetPitch

\phi

roll angle

GetRoll

r

body yaw rate

GetYawRate

q

body pitch rate

GetPitchRate

p

body roll rate

GetRollRate

\overline{c}

wing chord

wing_chord_ft

\overline{s}

wing span

wing_span_ft

t

time

\vec{x}

inertial position

\vec{v}

inertial velocity

\vec{a}

inertial acceleration

\vec{\omega}

rotational velocity

\vec{\alpha}

rotational acceleration

\hat{v}

velocity unit vector in inertial coordinates

\bold{R}

direction cosine matrix transforming a vector in the inertial frame to the body frame

m

mass

mass

GetEmptyWeight

\bold{I}

moment of inertia matrix

moment_of_inertia_ixx, moment_of_inertia_iyy, moment_of_inertia_izz

\vec{F}_{T,i}

total force in the inertial frame

\vec{M}_T

total moment

\vec{F}_{T,b}

total force in the body frame

\vec{F}_a

aerodynamic force in the body frame

\vec{F}_g

gravitational force in the body frame

\vec{F}_p

force from propulsion in the body frame

\vec{F}_l

force from the landing gear in the body frame

\vec{M}_a

moment produced from aerodynamics

\vec{r}_{cm/a}

center of mass relative to the aerodynamic reference point

center_of_mass_x1, center_of_mass_y1, center_of_mass_z1

GetCgX, GetCgY, GetCgZ

\vec{g}_i

gravitational vector in the inertial frame

\tilde{q}

dynamic pressure

A_r

reference area

reference_area

\gamma

area multiplier

C_L

lift coefficient

cL_alpha_beta_mach_table

C_{Lq}

lift coefficient contribution due to pitch rate, \frac{dC_L}{dq}

cLq_alpha_mach_table

C_{L\dot{\alpha}}

lift coefficient contribution due to the rate of change of angle of attack, \frac{dC_L}{d\dot{\alpha}}

cL_alphadot_alpha_mach_table

C_d

drag coefficient

cd_alpha_beta_mach_table

C_Y

side force coefficient

cy_alpha_beta_mach_table

C_{Yr}

side force coefficient contribution due to yaw rate, \frac{dC_Y}{dr}

cyr_beta_mach_table

C_{Y\dot{\beta}}

side force coefficient contribution due to the rate of change of side slip, \frac{dC_Y}{d\dot{\beta}}

cy_betadot_beta_mach_table

C_n

yaw moment coefficient

cn_alpha_beta_mach_table

C_{nr}

yaw moment coefficient contribution due to yaw rate, \frac{dC_n}{dr}

cnr_mach_table

C_{np}

yaw moment coefficient contribution due to roll rate, \frac{dC_n}{dp}

cnp_mach_table

C_{n\dot{\beta}}

yaw moment coefficient contribution due to the rate of change of side slip, \frac{dC_n}{d\dot{\beta}}

cn_betadot_mach_table

C_m

pitch moment coefficient

cm_alpha_beta_mach_table

C_{mq}

pitch moment coefficient contribution due to pitch rate, \frac{dC_m}{dq}

cmq_mach_table

C_{mp}

pitch moment coefficient contribution due to roll rate, \frac{dC_m}{dp}

cmp_mach_table

C_{m\dot{\alpha}}

pitch moment coefficient contribution due to the rate of change of angle of attack, \frac{dC_m}{d\dot{\alpha}}

cm_alphadot_mach_table

C_l

roll moment coefficient

cl_alpha_beta_mach_table

C_{lp}

roll moment coefficient contribution due to roll rate, \frac{dC_l}{dp}

clp_mach_table

C_{lr}

roll moment coefficient contribution due to yaw rate, \frac{dC_l}{dr}

clr_mach_table

C_{lq}

roll moment coefficient contribution due to pitch rate, \frac{dC_l}{dq}

clq_mach_table

C_{l\dot{\alpha}}

roll moment coefficient contribution due to the rate of change of angle of attack, \frac{dC_l}{d\dot{\alpha}}

cl_alphadot_mach_table

C_{l\dot{\beta}}

roll moment coefficient contribution due to the rate of change of side slip, \frac{dC_l}{d\dot{\beta}}

cl_betadot_mach_table

\hat{l}

lift unit vector in body coordinates

\hat{d}

drag unit vector in body coordinates

\hat{s}

side force unit vector in body coordinates

T

engine thrust, dependent on the engine type

\hat{t}

thrust unit vector in the body frame

Equations of Motion

The inertial acceleration is found by dividing the total inertial force by the total mass. In a similar manner, the rotational acceleration is found by dividing the total moment by the total moment of inertia. These accelerations, given by (1), are computed using the average force and moment between two time steps, as shown in equation (2).

(1)\vec{a} &= \frac{\overline{\vec{F}}_{T,i}}{m} \vec{\alpha} &= \bold{I}^{-1}\overline{\vec{M}}_T

(2)\overline{\vec{F}}_{T,i} = \frac{1}{2} [\vec{F}_{T,i}]_{t} + \frac{1}{2} [\vec{F}_{T,i}]_{t+\Delta t} \overline{\vec{M}}_{T} = \frac{1}{2} [\vec{M}_{T}]_{t} + \frac{1}{2} [\vec{M}_{T}]_{t+\Delta t}

It should be noted that one of the simplifications of P6DOF is that the moment of inertia from (1) is a diagonal matrix given by (3).

(3)\bold{I} = \left[ {\begin{array}{rrr} I_{xx} & 0 & 0 \\ 0 & I_{yy} & 0 \\ 0 & 0 & I_{zz} \end{array}} \right]

First, the force [\vec{F}_{T,i}]_{t} and moment [\vec{M}_{T}]_{t} are computed at the current time step, which are used to calculate the linear and angular acceleration. Then, the states are propagated forward in time using (4). The force and moment are then computed at this new state, producing [\vec{F}_{T,i}]_{t+\Delta t} and [\vec{M}_{T}]_{t+\Delta t}. The average force and the average moment are computed using (2). These are then used to compute the linear and angular accelerations in (1). The states are then updated using (4).

(4)\vec{x}_{t+\Delta t} &= \vec{x}_t + \vec{v}_t \Delta t + \frac{1}{2} \vec{a}_t {\Delta t}^2 \vec{v}_{t+\Delta t} &= \vec{v}_t + \vec{a}_t \Delta t \vec{Q}_{t+\Delta t} &= \vec{Q}_t + \vec{\dot{Q_t}} \Delta t \vec{\omega}_{t+\Delta t} &= \vec{\omega}_t + \vec{\alpha}_t \Delta t

The body angular rates \vec{\omega}, quaternion \vec{Q}, and rate quaternion \vec{\dot{Q}} are defined in (5).

(5)\vec{\omega} &= \left[ {\begin{array}{ccc} r & q & p \end{array}} \right]^\top \vec{Q} &= \left[ {\begin{array}{c} \cos{\frac{\psi}{2}} \cos{\frac{\theta}{2}} \cos{\frac{\phi}{2}} + \sin{\frac{\psi}{2}} \sin{\frac{\theta}{2}} \sin{\frac{\phi}{2}} \\ \cos{\frac{\psi}{2}} \cos{\frac{\theta}{2}} \sin{\frac{\phi}{2}} - \sin{\frac{\psi}{2}} \sin{\frac{\theta}{2}} \cos{\frac{\phi}{2}} \\ \cos{\frac{\psi}{2}} \sin{\frac{\theta}{2}} \cos{\frac{\phi}{2}} + \sin{\frac{\psi}{2}} \cos{\frac{\theta}{2}} \cos{\frac{\phi}{2}} \\ \sin{\frac{\psi}{2}} \cos{\frac{\theta}{2}} \cos{\frac{\phi}{2}} - \cos{\frac{\psi}{2}} \sin{\frac{\theta}{2}} \sin{\frac{\phi}{2}} \end{array}} \right] \vec{\dot{Q_t}} &= \frac{1}{2} \left[ {\begin{array}{cccc} 0 & -r & -q & -p \\ r & 0 & p & -q \\ q & -p & 0 & r \\ p & q & -r & 0 \end{array}} \right] \vec{Q}_t

If the quaternion \vec{Q} is defined by \left[ {\begin{array}{cccc}a & b & c & d\end{array}} \right]^\top, then the corresponding direction cosine matrix is computed using (6).

(6)\bold{R} = \left[ {\begin{array}{ccc} a^2 + b^2 - c^2 - d^2 & 2(bc + ad) & 2(bd - ac)\\ 2(bc - ad) & a^2 - b^2 + c^2 - d^2 & 2(cd + ab) \\ 2(bd + ac) & 2(cd - ab) & a^2 - b^2 - c^2 + d^2 \end{array}} \right]

Forces and Moments

The forces are calculated in the body frame, and then converted to the inertial frame to use in equation (1). This conversion is accomplished through equation (7).

(7)\vec{F}_{T,i} = \bold{R}^{-1}\vec{F}_{T,b}

where \bold{R} is defined as a 3-2-1 rotation matrix (yaw-pitch-roll) shown in equation (8).

(8)\bold{R} = \left[ {\begin{array}{rrr} \cos(\theta)\cos(\psi) & \cos(\theta)\sin(\psi) & -\sin(\theta) \\ \sin(\phi)\sin(\theta)\cos(\psi) - \cos(\phi)\sin(\psi) & \sin(\phi)\sin(\theta)\sin(\psi) + \cos(\phi)\cos(\psi) & \sin(\phi)\cos(\theta) \\ \cos(\phi)\sin(\theta)\cos(\psi) + \sin(\phi)\sin(\psi) & \cos(\phi)\sin(\theta)\sin(\psi) - \sin(\phi)\cos(\psi) & \cos(\phi)\cos(\theta) \\ \end{array}} \right]

The total body force, \vec{F}_{t,b} from equation (7) is calculated as the sum of the aerodynamic, gravitational, and propulsion forces in the body frame. It also includes, if applicable, the force from landing gear.

\vec{F}_{T,b} = \vec{F}_a + \vec{F}_g + \vec{F}_p + \vec{F}_l

The total moment is given by equation (9).

(9)\vec{M}_T = \vec{M}_a + \vec{r}_{cm/a} \times (\vec{F}_a + \vec{F}_p + \vec{F}_l)

The body frame gravitational force is calculated in (10).

(10)\vec{F}_g = m\bold{R}\vec{g}_i

The body frame propulsion force is given by (11).

(11)\vec{F}_p = T\hat{t}_b

The force and moment produced from aerodynamics is given by (12) and (13).

(12)\vec{F}_a = \tilde{q}A_r\gamma[(C_L(\alpha, \beta, M) + k_{Lq}C_{Lq}(\alpha, M) + k_{L\dot{\alpha}}C_{L\dot{\alpha}}(\alpha, M))&\hat{l} \\ + (C_Y(\alpha, \beta, M) + k_{Yr}C_{Yr}(\beta, M) + k_{Y\dot{\beta}}C_{Y\dot{\beta}}(\beta, M))&\hat{s} \\ + C_d(\alpha, \beta, M)&\hat{d}]

(13)\vec{M}_a = \tilde{q}A_r[ (C_l(\alpha, \beta, M) + k_{lp}C_{lp}(M) + k_{lr}C_{lr}(M) + k_{lq}C_{lq}(M) + k_{l\dot{\alpha}}C_{l\dot{\alpha}}(M) + k_{l\dot{\beta}}C_{l\dot{\beta}}(M))&\hat{i} \\ + (C_m(\alpha, \beta, M) + k_{mq}C_{mq}(M) + k_{mp}C_{mp}(M) + k_{m\dot{\alpha}}C_{m\dot{\alpha}}(M))&\hat{j} \\ + (C_n(\alpha, \beta, M) + k_{nr}C_{nr}(M) + k_{np}C_{np}(M) + k_{n\dot{\beta}}C_{n\dot{\beta}}(M))&\hat{k}]

The vectors \hat{i}, \hat{j}, \hat{k} in (13) are the unit vectors that form the standard orthonormal basis in the body frame.

The vectors \hat{l}, \hat{d}, \hat{s} in (12) are the lift, drag, and side force unit vectors in the body frame, and are computed using (14).

(14)\hat{l} &= \hat{j} \times \bold{R}\hat{v} \hat{d} &= - \bold{R}\hat{v} \hat{s} &= \hat{l} \times \hat{d}

Reduced Frequency

The reduced frequencies 1, denoted by k_{\{x\}}, in (12) and (13) are calculated using (15). If use_reduced_frequency is set to false, then the corresponding angular rates are used instead of the non-dimensional reduced frequencies, and the stability derivative tables should be adjusted accordingly.

(15)k_{Lq} &= \frac{\overline{c} q}{2 \lVert \vec{v} \rVert} k_{L\dot{\alpha}} &= \frac{\overline{c}\dot{\alpha}}{2 \lVert \vec{v} \rVert} k_{Yr} &= \frac{\overline{s} r}{2 \lVert \vec{v} \rVert} k_{Y\dot{\beta}} &= \frac{\overline{s} \dot{\beta}}{2 \lVert \vec{v} \rVert} k_{lp} &= \frac{\overline{s} p}{2 \lVert \vec{v} \rVert} k_{lr} &= \frac{\overline{s} r}{2 \lVert \vec{v} \rVert} k_{lq} &= \frac{\overline{s} q}{2 \lVert \vec{v} \rVert} k_{mq} &= \frac{\overline{c} q}{2 \lVert \vec{v} \rVert} k_{mp} &= \frac{\overline{c} p}{2 \lVert \vec{v} \rVert} k_{m\dot{\alpha}} &= \frac{\overline{c} \dot{\alpha}}{2 \lVert \vec{v} \rVert} k_{nr} &= \frac{\overline{s} r}{2 \lVert \vec{v} \rVert} k_{np} &= \frac{\overline{s} p}{2 \lVert \vec{v} \rVert} k_{n\dot{\beta}} &= \frac{\overline{s} \dot{\beta}}{2 \lVert \vec{v} \rVert}

The terms \overline{c} and \overline{s} refer to the wing chord and span, respectively.

References

1

Jenkins, J. “Dynamic Stability Derivatives”, Air Force Research Laboratory - Aerodynamic Technology Branch Aerospace Vehicles Division, June 2015. AFRL-RQ-WP-TR-2015-0141.