2.1 Reference frames and kinematic transformations

The derivation of a control-oriented nonlinear model for the F-35B aircraft necessitates a rigorous definition of the coordinate systems and the mathematical transformations that link them. The inertial reference frame and the body-fixed reference frame are the two main reference frames used in this study to characterise the aircraft’s motion.

When selecting a coordinate system for aircraft or airborne systems in general, there are a few alternatives. To track the absolute trajectory of the aircraft, particularly during the critical transition to landing on an inclined surface, the North-East-Down (NED) system is adopted as the inertial frame [Reference Koks25]. The position of the aircraft’s centre of gravity (CG) in this frame is denoted by the vector:

(1) \begin{align}\boldsymbol{r}={\left[\begin{array}{ccc}X \quad Y \quad Z\end{array}\right]}^{\mathrm{T}}\end{align}

where $Z$ is positive downward according to standard aerospace conventions. The body-axis convention is shown in Fig. 1, with ${x}_{b}$ forward, ${y}_{b}$ right and ${z}_{b}$ downward. The state of the aircraft within this local frame is defined by its linear velocities and angular rates are defined in Equations (2) and (3).

(2) \begin{align}{\mathbf{v}}_{b}={\left[\begin{array}{ccc}u \quad v {}\quad w\end{array}\right]}^{\mathrm{T}} \end{align}

(3) \begin{align}\begin{array}{c}{\boldsymbol{\omega }}_{b}={\left[\begin{array}{ccc}p \quad q \quad r\end{array}\right]}^{\mathrm{T}} \end{array}\end{align}

The orientation, or attitude, of the aircraft relative to the NED frame is parameterised by the Euler angles:

(4) \begin{align}\begin{array}{c}\boldsymbol{\eta} ={\left[\begin{array}{ccc}\phi \quad \theta \quad \psi \end{array}\right]}^{\mathrm{T}}\end{array}\end{align}

representing the roll, pitch and yaw, respectively. The mapping of translational velocities from the body frame to the inertial frame is achieved through a direction cosine matrix (DCM), denoted as ${\boldsymbol{R}}_{b}^{i}\left(\phi ,\theta ,\psi \right)$ . Following the standard ZYX (3-2-1) Euler angle sequence, the rotation matrix is constructed by performing successive rotations about the yaw, pitch and roll axes. This matrix ensures that any vector in the body frame ( ${\boldsymbol{a}}_{\boldsymbol{b}}$ ) can be projected into the inertial frame ( ${\boldsymbol{a}}_{\boldsymbol{i}}$ ) via the given relation:

(5) \begin{align}\begin{array}{c}{\boldsymbol{a}}_{\boldsymbol{i}}={\boldsymbol{R}}_{b}^{i}{\boldsymbol{a}}_{\boldsymbol{b}}\end{array}\end{align}

The explicit form of the DCM used in the simulation framework is given by

(6) \begin{align}\begin{array}{c}{\boldsymbol{R}}_{b}^{i}=\left[\begin{array}{c@{\quad}c@{\quad}c}\mathrm{cos}\psi \mathrm{cos}\theta & {} \mathrm{cos}\psi \mathrm{sin}\theta \mathrm{sin}\psi -\mathrm{cos}\phi \mathrm{sin}\psi & \mathrm{cos}\psi \mathrm{sin}\theta \mathrm{cos}\phi +\mathrm{sin}\psi \mathrm{sin}\phi \\ \mathrm{cos}\theta \mathrm{sin}\psi & \mathrm{sin}\psi \mathrm{sin}\theta \mathrm{sin}\phi +\mathrm{cos}\psi \mathrm{cos}\phi & \mathrm{sin}\psi \mathrm{sin}\theta \mathrm{cos}\phi -\mathrm{cos}\psi \mathrm{sin}\phi \\ -\mathrm{sin}\theta & \mathrm{cos}\theta \mathrm{sin}\phi & \mathrm{cos}\theta \mathrm{cos}\phi \end{array}\right]\end{array}\end{align}

Since the transformation matrix is orthonormal, the inverse transformation (inertial to body) is simply its transpose: ${\boldsymbol{R}}_{i}^{b}={\left({\boldsymbol{R}}_{b}^{i}\right)}^{\mathrm{T}}$ . The translational motion of the CG is described in the inertial frame. The time derivatives of the inertial position, which define the aircraft’s trajectory, are obtained by multiplying the body velocity vector ${\mathbf{v}}_{b}$ by the DCM:

(7) \begin{align}\begin{array}{c}\dot{\boldsymbol{r}}=\left[\begin{array}{c}\,\dot{\!X}\\ \,\dot{\!Y}\\ \,\dot{\!Z}\end{array}\right]={\boldsymbol{R}}_{b}^{i}\left[\begin{array}{c}u\\ v\\ w\end{array}\right]\end{array}\end{align}

Expanding this matrix-vector product yields the explicit scalar kinematic equations used in the computational model:

(8) \begin{align}\begin{array}{c}\,\dot{\!X}=u\left(\mathrm{cos}\theta \mathrm{cos}\psi \right)+v\left(\mathrm{sin}\phi \mathrm{sin}\theta \mathrm{cos}\psi -\mathrm{cos}\phi \mathrm{sin}\psi \right)\\ +\, w\left(\mathrm{cos}\phi \mathrm{sin}\theta \mathrm{cos}\psi +\mathrm{sin}\phi \mathrm{sin}\psi \right)\end{array}\end{align}

(9) \begin{align}\begin{array}{c}\,\dot{\!Y}=u\left(\mathrm{cos}\mathrm{sin}\psi \right)+v\left(\mathrm{sin}\phi \mathrm{sin}\theta \mathrm{sin}\psi +\mathrm{cos}\phi \mathrm{cos}\psi \right)\\ +\, w\left(\mathrm{cos}\mathrm{sin}\theta \mathrm{sin}\psi -\mathrm{sin}\phi \mathrm{cos}\psi \right)\end{array}\end{align}

(10) \begin{align}\begin{array}{c}\,\dot{\!Z}=-u\mathrm{sin}\theta +v\mathrm{sin}\phi \mathrm{cos}\theta +w\mathrm{cos}\phi \mathrm{cos}\theta \end{array}\end{align}

A common misconception in simplified flight models is the assumption that body angular rates ( $p$ , $q$ , $r$ ) are identical to the rates of change of the Euler angles ( $\dot{\phi }$ , $\dot{\theta }$ , $\dot{\psi }$ ). In a rotating 3D environment, these rates are related through a nonlinear transformation. This relationship is governed by the transformation matrix $\boldsymbol{T}(\phi ,\theta )$ :

(11) \begin{align}\begin{array}{c}\dot{\boldsymbol{\eta }}=\boldsymbol{T}(\phi ,\theta ){\boldsymbol{\omega }}_{b}\end{array}\end{align}

where the expanded form is defined as:

(12) \begin{align}\begin{array}{c}\left[\begin{array}{c}\dot{\phi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right]=\left[\begin{array}{c@{\quad}c@{\quad}c}1 & \mathrm{sin}\phi \mathrm{t}\mathrm{a}\mathrm{n}\theta & \mathrm{cos}\phi \mathrm{t}\mathrm{a}\mathrm{n}\theta \\ 0 & \mathrm{cos}\phi & -\mathrm{sin}\phi \\ 0 {}& \mathrm{sin}\phi \mathrm{s}\mathrm{e}\mathrm{c}\theta {}& \mathrm{cos}\phi \mathrm{s}\mathrm{e}\mathrm{c}\theta \end{array}\right]\left[\begin{array}{c}p\\ q\\ r\end{array}\right]\end{array}\end{align}

This relationship leads to the following set of first-order differential equations used to propagate the aircraft’s attitude:

(13) \begin{align}\begin{array}{c}\dot{\phi }=p+tan\theta \left(q\mathrm{sin}\phi +r\mathrm{cos}\phi \right)\end{array}\end{align}

(14) \begin{align}\begin{array}{c}\dot{\theta }=qcos\phi -rsin\phi \end{array}\end{align}

(15) \begin{align}\dot{\psi }=\dfrac{q\mathrm{sin}\phi +r\mathrm{cos}\phi }{\mathrm{cos}\theta }\end{align}

The use of Euler angles introduces a mathematical singularity known as gimbal lock when the pitch angle ${\theta }$ reaches ±90^∘, causing the $\mathrm{s}\mathrm{e}\mathrm{c}\theta $ and $\mathrm{t}\mathrm{a}\mathrm{n}\theta $ terms to approach infinity. While the typical landing envelopes for the F-35B do not involve such extreme pitch attitudes, the simulation architecture incorporates a numerical guard to maintain stability. Specifically, if the value of $\left|\mathrm{cos}\theta \right|$ falls below a threshold of ${10}^{-3}$ , it is constrained to $\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\mathrm{cos}\theta \right)\cdot {10}^{-3}$ , preventing division-by-zero errors and ensuring the robustness of the trajectory propagation.

In addition, the angular velocity vector in the body frame, ${\boldsymbol{\omega }}_{b}={\left[\begin{array}{ccc}p \quad q\quad r\end{array}\right]}^{\mathrm{T}}$ , is related to the Euler angle rates $\dot{\boldsymbol{\eta }}$ via the kinematic transformation matrix $\boldsymbol{W}\left(\phi ,\theta \right)$ :

(16) \begin{align}\begin{array}{c}{\boldsymbol{\omega }}_{b}=\boldsymbol{W}\left(\phi ,\theta \right)\dot{\boldsymbol{\eta }}\end{array}\end{align}

where:

(17) \begin{align}\begin{array}{c}\boldsymbol{W}\left(\phi ,\theta \right)=\left[\begin{array}{c@{\quad}c@{\quad}c}1 & 0 & -\mathrm{sin}\theta \\ 0 {}& \mathrm{cos}\phi & \mathrm{sin}\phi \mathrm{cos}\theta \\ 0 & -\mathrm{sin}\phi & \mathrm{cos}\phi \mathrm{cos}\theta \end{array}\right]\end{array}\end{align}

Equation (16) implies the familiar scalar relations:

(18) \begin{align} p=\dot{\phi }-\dot{\psi }\mathrm{sin}\theta \end{align}

(19) \begin{align}\begin{array}{c}q=\dot{\theta }\mathrm{cos}\phi +\dot{\psi }\mathrm{sin}\phi \mathrm{cos}\theta \end{array}\end{align}

(20) \begin{align}\begin{array}{c}r=-\dot{\theta }\mathrm{sin}\phi +\dot{\psi }\mathrm{cos}\phi \mathrm{cos}\theta \end{array}\end{align}

These expressions highlight how the body rates are synthesised from the individual Euler rate components, adjusted by the current orientation of the aircraft. Beyond simple kinematic mapping, the $\boldsymbol{W}\left(\phi ,\theta \right)$ matrix is fundamental for calculating the rotational kinetic energy ( ${T}_{rot}$ ) within the Lagrangian framework. By substituting the relationship ${\boldsymbol{\omega }}_{b}=\boldsymbol{W}\left(\phi ,\theta \right)\dot{\boldsymbol{\eta }}$ into the standard rotational energy equation, the energy is expressed purely in terms of Euler rates:

(21) \begin{align}\begin{array}{c}{T}_{rot}=\dfrac{1}{2}{\boldsymbol{\omega }}_{b}^{\mathrm{T}}\boldsymbol{I}{\boldsymbol{\omega }}_{b}=\dfrac{1}{2}{\dot{\boldsymbol{\eta }}}^{\mathrm{T}}\left({\boldsymbol{W}}^{\mathrm{T}}\left(\phi ,\theta \right)\mathrm{}\boldsymbol{I}\mathrm{}\boldsymbol{W}\left(\phi ,\theta \right)\right)\dot{\boldsymbol{\eta }}\end{array}\end{align}

This leads to the definition of the Euler-rate inertia matrix, $\boldsymbol{J}\left(\phi ,\theta \right)={\boldsymbol{W}}^{\mathrm{T}}\left(\phi ,\theta \right)\mathrm{}\boldsymbol{I}\mathrm{}\boldsymbol{W}\left(\phi ,\theta \right)$ which retains all cross-inertia couplings and allows for the derivation of the equations of motion using generalised coordinates. This dual-representation capability (using both $\boldsymbol{T}\left(\phi ,\theta \right)$ and $\boldsymbol{W}\left(\phi ,\theta \right)$ matrices) ensures numerical stability and theoretical consistency across the different modules of the F-35B simulation.

2.2 Governing equations of motion

The dynamics are resolved by aggregating all generalised forces and moments acting on the CG. These include propulsive thrust ${\boldsymbol{F}}_{th,i}$ , ${\boldsymbol{M}}_{th,i}$ , gravitational effects ${\boldsymbol{F}}_{g,i}$ , aerodynamic loads ${\boldsymbol{F}}_{aero,i}$ , ${\boldsymbol{M}}_{aero,i}$ , landing gear reactions ${\boldsymbol{F}}_{gear,i}$ , ${\boldsymbol{M}}_{gear,i}$ and external disturbances ${\boldsymbol{F}}_{dist,i}$ , ${\boldsymbol{M}}_{dist,i}$ where $i=x,y,z$ . In the nominal simulations, the external disturbance force and moment terms are set to zero. In the additional disturbance case presented in Section 5.2, bounded gust-equivalent disturbance forces are introduced through the same external-force channels.

Newton’s Second Law governs translational motion. However, because the body frame is a non-inertial rotating reference, the acceleration must account for the Coriolis (transport) acceleration the ${\boldsymbol{\omega }}_{b}\times {\mathbf{v}}_{b}$ term. The fundamental equation is expressed as:

(22) \begin{align}\begin{array}{c}m({\dot{\mathbf{v}}}_{b}+{\boldsymbol{\omega }}_{b}\times {\mathbf{v}}_{b})={\boldsymbol{F}}_{total}\end{array}\end{align}

To determine the state derivatives ( $\dot{u}$ , $\dot{v}$ , $\dot{w}$ ), we first synthesise the total body-frame forces. The full form of the ${\boldsymbol{F}}_{total}$ vector is as follows:

(23) \begin{align}\begin{array}{c}{\boldsymbol{F}}_{total}=\left[\begin{array}{c}{F}_{total,x}\\ {F}_{total,y}\\ {F}_{total,z}\end{array}\right]\end{array}\end{align}

The F-35B propulsion system is unique, utilising a lift fan ( ${F}_{\,f}$ ), a rear swivel nozzle ( ${F}_{r}$ ) and two wing roll posts ( ${F}_{w1}$ , ${F}_{w2}$ ). The thrust directions for the lift fan and rear nozzle are variable, controlled by the deflection angles ${\delta }_{\,f}$ and ${\delta }_{r}$ , respectively. In the hover configuration ( $\delta =0$ ), the thrust is directed vertically along the $+{z}_{b}$ axis. Conversely, for $\delta \gt 0$ , the thrust vector tilts towards the $+{x}_{b}$ axis, thereby generating a forward acceleration component alongside the vertical lift. The thrust from the front and rear nozzles is decomposed based on their deflection angles ( ${\delta }_{\,f}$ , ${\delta }_{r}$ ):

(24) \begin{align}\begin{array}{c}{F}_{th,x}={F}_{\,f}\,\ \mathrm{sin}\ \,{\delta }_{\,f}+{F}_{r}\ \,\mathrm{sin}\,\ {\delta }_{r}\end{array}\end{align}

(25) \begin{align}\begin{array}{c}{F}_{th,y}=0\end{array}\end{align}

(26) \begin{align}\begin{array}{c}{{F}_{th,z}=F}_{\,f}\,\ \mathrm{cos}\,\ {\delta}_{\,f}+{F}_{r}\, \ \mathrm{cos}\, \ {\delta }_{r}+{F}_{w1}+{F}_{w2}\end{array}\end{align}

Note that the lateral thrust component ${F}_{th,y}$ is considered zero in this configuration.

Gravity acts along the $+Z$ axis in the NED frame. To map it into the body frame, we use the Euler angles $\phi $ and $\theta $ :

(27) \begin{align}\begin{array}{c}{F}_{g,x}=-mgsin\theta \end{array}\end{align}

(28) \begin{align}\begin{array}{c}{F}_{g,y}=mgsin\phi cos\theta \end{array}\end{align}

(29) \begin{align}\begin{array}{c}{F}_{g,z}=mgcos\phi cos\theta \end{array}\end{align}

Using ${\boldsymbol{F}}_{aero}$ and ${\boldsymbol{F}}_{gear}$ , which will be introduced later in the article, the total force expression can be obtained as follows.

(30) \begin{align}\begin{array}{c}{F}_{total,i}={F}_{th,i}+{F}_{g,i}+{F}_{aero,i}+{F}_{gear,i}+{F}_{dist,i}\ \text{where}\ i\in \left\{x,y,z\right\}\end{array}\end{align}

By expanding the cross-product ${\boldsymbol{\omega }}_{b}\times {\mathbf{v}}_{b}$ we arrive at the explicit differential equations for the linear velocities:

(31) \begin{align}\begin{array}{c}\dot{u}=\dfrac{1}{m}{F}_{total,x}-\left(qw-rv\right)\end{array}\end{align}

(32) \begin{align}\begin{array}{c}\dot{v}=\dfrac{1}{m}{F}_{total,y}-\left(ru-pw\right)\end{array}\end{align}

(33) \begin{align}\begin{array}{c}\dot{w}=\dfrac{1}{m}{F}_{total,z}-\left(pv-qu\right)\end{array}\end{align}

The rotational motion is described by Euler’s equations, which relate the total external moment ${\boldsymbol{M}}_{total}$ to the rate of change of angular momentum. ${\boldsymbol{M}}_{total}$ is the sum of moments from thrust, aerodynamics, the landing gear and disturbances. The full form of the ${\boldsymbol{M}_{total}}$ vector is as follows:

(34) \begin{align}\begin{array}{*{20}{c}} {\boldsymbol{M}_{total}} = \left[ {\begin{array}{*{20}{c}}{M_{total,x}}\\{{M_{total,y}}}\\{{M_{total,z}}}\end{array}} \right] \end{array}\end{align}

The rigid-body inertia tensor $\boldsymbol{I}$ includes full cross-coupling terms to ensure fidelity:

(35) \begin{align}\begin{array}{*{20}{c}}\boldsymbol{I} = \left[ {\begin{array}{c@{\quad}c@{\quad}c}{{I_{xx}}} & { - {I_{xy}}} {}&{ - {I_{xz}}}\\{ - {I_{xy}}}& {{I_{yy}}} & { - {I_{yz}}}\\{ - {I_{xz}}} & { - {I_{yz}}}& {{I_{zz}}}\end{array}} \right] \end{array}\end{align}

In this matrix, the diagonal elements ( ${I_{xx}},{I_{yy}},{I_{zz}}$ ) represent the moments of inertia, quantifying the aircraft’s resistance to rotational acceleration about the $x$ , $y$ and $z$ axes, respectively. The off-diagonal terms (e.g. ${I_{xy}}$ , ${I_{xz}}$ ) correspond to the products of inertia, which describe the mass symmetry relative to the orthogonal planes. Due to the reciprocal property of mass distribution, the tensor is symmetric, implying that ${I_{xy}} = {I_{yx}}$ , ${I_{xz}} = {I_{zx}}$ , and ${I_{yz}} = {I_{zy}}$ . For an aircraft with geometric symmetry across the $xz$ -plane, the products of inertia ${I_{xy}}$ and ${I_{yz}}$ are typically negligible and assumed to be zero, leaving ${I_{xz}}$ as the primary coupling term.

The governing rotational equation for rotational dynamics with inertia coupling is:

(36) \begin{align}\begin{array}{*{20}{c}}\boldsymbol{I}{{\dot {\boldsymbol{\omega}} }_b} + {\boldsymbol{\omega} _b} \times \left( {\boldsymbol{I}{\boldsymbol{\omega} _b}} \right) = {\boldsymbol{M}_{total}} \end{array}\end{align}

The total thrust and moments are generated by four distinct effectors: the lift fan ( ${F_f}$ ), the rear swivel nozzle ( ${F_r}$ ) and two wing roll posts ( ${F_{w1}}$ , ${F_{w2}}$ ). Moments are generated by the cross-product of the displacement vectors ${\boldsymbol{r}_k}$ (from CG to each nozzle) and their respective thrust vectors ${\boldsymbol{F}_k}$ (where $k \in \{f,r, w1, w2\}$ ). All displacement vectors are defined in Equations (37)–(40):

(37) \begin{align}\begin{array}{*{20}{c}}{\boldsymbol{r}_f} = \left[ {\begin{array}{*{20}{c}}{{L_f}}\\0\\0\end{array}} \right] \end{array}\end{align}

(38) \begin{align}\begin{array}{*{20}{c}}{\boldsymbol{r}_r} = \left[ {\begin{array}{*{20}{c}}{ - {L_r}}\\0\\0\end{array}} \right] \end{array}\end{align}

(39) \begin{align}\begin{array}{*{20}{c}}{\boldsymbol{r}_{w1}} = \left[ {\begin{array}{*{20}{c}}{{L_{w}}}\\{{L_b}}\\0\end{array}} \right] \end{array}\end{align}

(40) \begin{align}\begin{array}{*{20}{c}}{\boldsymbol{r}_{w2}} = \left[ {\begin{array}{*{20}{c}}{ - {L_w}}\\{ - {L_b}}\\0\end{array}} \right] \end{array}\end{align}

${L_f}$ , ${L_r}$ , ${L_w}$ , and ${L_b}$ values are given in Table 1 and are represented schematically in Fig. 1.

Table 1.

Aircraft geometry, mass and environment parameters

Table 1. Long description

A table with 14 rows and 4 columns detailing various parameters related to aircraft geometry, mass, and environment. The columns are labeled Symbol/variable, Description, Recommended value, and Unit. The table includes parameters such as gravitational acceleration, sea level air density, aircraft mass, wing reference area, wingspan, mean aerodynamic chord, roll-axis inertia, pitch-axis inertia, yaw-axis inertia, cross inertia term, front nozzle x-offset from CG, rear nozzle x-offset from CG, roll-post x-offset from CG, and roll-post y-offset from CG. Each row provides the symbol or variable, its description, the recommended value, and the unit of measurement. For example, the gravitational acceleration is listed as 9.81 m/s2, sea level air density as 1.225 kg/m3, and aircraft mass as 20,000 kg. The table provides a comprehensive overview of the key parameters involved in aircraft design and performance.

Although the geometric offset parameters are denoted by (L)-based symbols for consistency with the original model notation, their coordinate directions are explicitly specified through the body-frame position vectors. Positive and negative signs in these vectors indicate whether the corresponding offset is measured along the positive or negative ( ${x_b}$ ), ( ${y_b}$ ) or ( ${z_b}$ ) direction. Table 1 also summarises the physical parameters of the proposed aircraft system. The thrust moment about the CG, which provides primary control over the roll and pitch channels, is computed as a sum of cross products:

(41) \begin{align}\begin{array}{*{20}{c}}{\boldsymbol{M}_{th}} = \mathop \sum \limits_{k \in \left\{ {\,f,r,w1,w2} \right\}} {\boldsymbol{r}_k} \times {\boldsymbol{F}_k} \end{array}\end{align}

Using Equations (37)–(41), the resulting moment components are:

(42) \begin{align}\begin{array}{*{20}{c}}{M_{th,x}} = {L_b}\left( {{F_{w1}} - {F_{w2}}} \right) \end{array}\end{align}

(43) \begin{align}\begin{array}{*{20}{c}}{M_{th,y}} = - {L_f}\left( {{F_f}{\rm{cos}}{\delta _f}} \right) + {L_r}\left( {{F_r}{\rm{cos}}{\delta _r}} \right) + {L_w}\left( {{F_{w1}} + {F_{w2}}} \right) \end{array}\end{align}

(44) \begin{align}\begin{array}{*{20}{c}}{M_{th,z}} \approx 0 \end{array}\end{align}

where yaw control is neglected at this stage.

The gyroscopic (Coriolis) coupling term $\boldsymbol{C} = {\boldsymbol{\omega} _b} \times \left( {\boldsymbol{I}{\boldsymbol{\omega} _b}} \right)$ represents the inertial moments arising from the rotation of the aircraft. For the F-35B, these terms capture how high-rate yaw or pitch can ‘spill over’ into other axes, a critical factor during aggressive landing corrections.

To solve for the angular accelerations ( $\dot p,\dot q,\dot r$ ), Equation (36) is rearranged, and the inertia matrix is inverted:

(45) \begin{align}\begin{array}{*{20}{c}}\boldsymbol{\dot{\omega}} = {\boldsymbol{I}^{ - 1}}\left( {\boldsymbol{M}_{total} - \boldsymbol{C}} \right) \end{array}\end{align}

For numerical robustness, the simulation framework includes a safety check on the determinant of $\boldsymbol{I}$ to prevent singularities, although the tensor is physically required to be positive definite.

The term $\boldsymbol{C} = {\boldsymbol{\omega} _b} \times \left( {\boldsymbol{I}{\boldsymbol{\omega} _b}} \right)$ represents the inertial moments arising from the rotation of the body-fixed frame. To avoid shorthand notation, we substitute the full components of the angular momentum directly. The resulting gyroscopic moment vector $\boldsymbol{C} = {\left[ {\begin{array}{*{20}{c}}{{C_x}} {}{{C_y}} {}{{C_z}}\end{array}} \right]^T}$ is derived as follows:

(46) \begin{align}\begin{array}{*{20}{c}}{C_x} = q\left( { - {I_{xz}}p - {I_{yz}}q + {I_{zz}}r} \right) - r\left( { - {I_{xy}}p + {I_{yy}}q - {I_{yz}}r} \right) \end{array}\end{align}

(47) \begin{align}\begin{array}{*{20}{c}}{C_y} = r\left( {{I_{xx}}p - {I_{xy}}q - {I_{xz}}r} \right) - p\left( { - {I_{xz}}p - {I_{yz}}q + {I_{zz}}r} \right) \end{array}\end{align}

(48) \begin{align}\begin{array}{*{20}{c}}{C_z} = p\left( { - {I_{xy}}p + {I_{yy}}q - {I_{yz}}r} \right) - q\left( {{I_{xx}}p - {I_{xy}}q - {I_{xz}}r} \right) \end{array}\end{align}

Here ${C_x}$ , ${C_y}$ and ${C_z}$ represent the roll axis, pitch axis and yaw axis gyroscopic term, respectively. These terms capture the spill-over effects where, for example, a simultaneous pitch ( $q$ ) and yaw ( $r$ ) rate generates a secondary rolling moment due to the distribution of the aircraft’s mass. By rearranging the fundamental equation, we isolate the product of the inertia tensor and the angular acceleration vector.

(49) \begin{align}\begin{array}{*{20}{c}}\left[ {\begin{array}{c@{\quad}c@{\quad}c}{{I_{xx}}} & { - {I_{xy}}} & { - {I_{xz}}}\\{ - {I_{xy}}} & {{I_{yy}}} & { - {I_{yz}}}\\{ - {I_{xz}}} & { - {I_{yz}}} & {{I_{zz}}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{\dot p}\\{\dot q}\\{\dot r}\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{M_{total,x}} - {C_x}}\\{{M_{total,y}} - {C_y}}\\{{M_{total,z}} - {C_z}}\end{array}} \right] \end{array}\end{align}

The instantaneous angular accelerations are obtained by multiplying the net effective moment by the inverse of the inertia tensor:

(50) \begin{align}\left[ {\begin{array}{*{20}{c}}{\dot p}\\{\dot q}\\{\dot r}\end{array}} \right] = {\left[ {\begin{array}{c@{\quad}c@{\quad}c}{{I_{xx}}}& { - {I_{xy}}}& { - {I_{xz}}}\\{ - {I_{xy}}}& {{I_{yy}}}& { - {I_{yz}}}\\{ - {I_{xz}}}& { - {I_{yz}}} & {{I_{zz}}}\end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}}{{M_{total,x}} - \left( { - q{I_{xz}}p - q{I_{yz}}q + q{I_{zz}}r + r{I_{xy}}p - r{I_{yy}}q + r{I_{yz}}r} \right)}\\{{M_{total,y}} - \left( {r{I_{xx}}p - r{I_{xy}}q - r{I_{xz}}r + p{I_{xz}}p + p{I_{yz}}q - p{I_{zz}}r} \right)}\\{{M_{total,z}} - \left( { - p{I_{xy}}p + p{I_{yy}}q - p{I_{yz}}r - q{I_{xx}}p + q{I_{xy}}q + q{I_{xz}}r} \right)}\end{array}} \right]\end{align}

For the F-35B, while ${I_{xy}}$ and ${I_{yz}}$ are often negligible due to lateral symmetry, the ${I_{xz}}$ term is retained as it provides critical coupling between the roll and yaw channels during asymmetric thrust conditions or crosswind landing manoeuvres. This complete 6-DOF dynamic model forms the basis for the subsequent development of aerodynamic and landing gear interaction modules.

2.3 Aerodynamic force and moment model

After the rigid-body equations of motion were derived, this section details the formulation of the aerodynamic forces and moments. For the F-35B, modelling these effects is particularly challenging during the STOVL phase, where the aircraft transitions from wing-borne flight to jet-borne hover. To capture these effects in a control-oriented manner, we employ a low-order aerodynamic surrogate model in which the aerodynamic coefficients are expressed as linear functions and normalised angular rates. This choice provides a practical balance between computational efficiency and physical interpretability, while remaining suitable for integration within the nonlinear 6-DOF simulation and subsequent controller design. The aerodynamic loads acting on the airframe are synthesised based on the aircraft’s interaction with the surrounding atmosphere. These loads are treated as non-conservative generalised forces ${\boldsymbol{F}_{aero}}$ and moments ${\boldsymbol{M}_{aero}}$ applied at the centre of gravity.

The aerodynamic model begins by defining the state of the airflow relative to the airframe. In the nominal simulations, the wind velocity is assumed to be zero; therefore, the body-fixed translational velocity components coincide with the air-relative velocity components. More generally, if an explicit wind field is introduced, the body-frame air-relative velocity vector should be computed as:

(51) \begin{align}\begin{array}{*{20}{c}}\boldsymbol{V}_{rel}^b = \left[ {\begin{array}{*{20}{c}}{{u_r}}\\{{v_r}}\\{{w_r}}\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}u\\v\\w\end{array}} \right] - {{\left( {\boldsymbol{R}_b^i} \right)}^{\mathrm{T}}}\boldsymbol{V}_w^n \end{array}\end{align}

where $\boldsymbol{V}_w^n$ is the wind velocity vector expressed in the NED inertial frame and ${\left( {\boldsymbol{R}_b^i} \right)^{\mathrm{T}}}$ maps this vector from the NED frame to the body frame. In the present nominal case, $\boldsymbol{V}_w^n = 0$ , so that ${u_r} = u$ , ${v_{}} = v$ and ${w_r} = w$ . The total airspeed is then defined from the air-relative velocity components as

(52) \begin{align}\begin{array}{*{20}{c}}{V_{air}} = \sqrt {u_r^2 + v_r^2 + w_r^2} \end{array}\end{align}

The orientation of the velocity vector relative to the fuselage is defined by the angle-of-attack ( $\alpha $ ) and the sideslip angle ( $\beta $ ):

(53) \begin{align}\begin{array}{*{20}{c}}\alpha = {\rm{atan}}2\left( {{w_r},{u_r}} \right) \end{array}\end{align}

(54) \begin{align}\begin{array}{*{20}{c}}\beta = {\rm{atan}}2\left( {{v_r},\sqrt {u_r^2 + w_r^2} } \right) \end{array}\end{align}

Finally, the dynamic pressure ( ${q_\infty }$ ), which scales all aerodynamic forces, is calculated as:

(55) \begin{align}\begin{array}{*{20}{c}}{q_\infty } = \dfrac{1}{2}\rho V_{air}^2 \end{array}\end{align}

This formulation explicitly distinguishes the aerodynamic relative velocity from the aircraft body-frame translational velocity. The no-wind nominal simulations reduce to the original expressions, whereas an external wind field can be incorporated through $\boldsymbol{V}_w^n$ . In the present study, the additional wind-gust case in Section 5.2 is introduced as a bounded gust-equivalent external disturbance force rather than as a full stochastic air-relative aerodynamic model.

The aerodynamic forces are first evaluated in the stability (wind) axes. The lift ( ${C_L}$ ), drag ( ${C_D}$ ) and side-force ( ${C_Y}$ ) coefficients follow a standard linear dependence on ${\rm{\alpha }}$ and sideslip, combined with a quadratic drag polar:

(56) \begin{align}\begin{array}{*{20}{c}}{C_L} = {C_{L0}} + {C_{L\alpha }}\alpha \end{array}\end{align}

(57) \begin{align}\begin{array}{*{20}{c}}{C_D} = {C_{D0}} + {K_{CD}}C_L^2 \end{array}\end{align}

(58) \begin{align}\begin{array}{*{20}{c}}{C_Y} = {C_{Y\beta }}\beta \end{array}\end{align}

where ${K_{CD}}$ represents the induced drag factor. The resulting stability-axis forces are calculated using the dynamic pressure ${q_\infty }$ and reference area $S$ :

(59) \begin{align}\begin{array}{*{20}{c}}{X_s} = - {q_\infty }S{C_D} \end{array}\end{align}

(60) \begin{align}\begin{array}{*{20}{c}}{Y_s} = {q_\infty }S{C_Y} \end{array}\end{align}

(61) \begin{align}\begin{array}{*{20}{c}}{Z_s} = - {q_\infty }S{C_L} \end{array}\end{align}

The body-fixed frame is then created by rotating these forces around the ${y_b}$ -axis. The body-frame aerodynamic force vector ${\boldsymbol{F}_{aero}} = {\left[ {\begin{array}{*{20}{c}}{{F_{aero,x}}}\quad {}{{F_{aero,y}}}\quad {{F_{aero,z}}}\end{array}} \right]^{\mathrm{T}}}$ is obtained as:

(62) \begin{align}\begin{array}{*{20}{c}}{F_{aero,x}} = {X_s}cos\alpha - {Z_s}sin\alpha \end{array}\end{align}

(63) \begin{align}\begin{array}{*{20}{c}}{F_{aero,y}} = {Y_s} \end{array}\end{align}

(64) \begin{align}\begin{array}{*{20}{c}}{F_{aero,z}} = {X_s}sin\alpha + {Z_s}cos\alpha \end{array}\end{align}

The aerodynamic moments are characterised by non-dimensional stability derivatives relative to the reference span ( $b$ ) and mean aerodynamic chord ( $c$ ). The formulation encompasses static stability terms and dynamic damping derivatives (dependent on normalised rates $p$ , $q$ , $r$ ):

(65) \begin{align} C_l = C_{l \beta}\beta + C_{lp} \frac{b}{2V_{air}}P + C_{lr} \frac{b}{2V_{air}}r \end{align}

(66) \begin{align}\begin{array}{*{20}{c}}{C_m} = {C_{m0}} + {C_{m\alpha }}\alpha + {C_{mq}}\frac{c}{{2{V_{air}}}}q \end{array}\end{align}

(67) \begin{align}\begin{array}{*{20}{c}}{C_n} = {C_{n{\rm{\beta }}}}\beta + {C_{np}}\frac{b}{{2{V_{air}}}}p + {C_{nr}}\frac{b}{{2{V_{air}}}}r \end{array}\end{align}

Finally, the body-frame aerodynamic moments ${\boldsymbol{M}_{aero}} = {\left[ {\begin{array}{*{20}{c}}{{M_{aero,x}}}\,\, {}{{M_{aero,y}}}\,\, {}{{M_{aero,z}}}\end{array}} \right]^{\mathrm{T}}}$ are scaled as follows:

(68) \begin{align}\begin{array}{*{20}{c}}{M_{aero,x}} = {q_\infty }Sb{C_l} \end{array}\end{align}

(69) \begin{align}\begin{array}{*{20}{c}}{M_{aero,y}} = {q_\infty }Sc{C_m} \end{array}\end{align}

(70) \begin{align}\begin{array}{*{20}{c}}{M_{aero,z}} = {q_\infty }Sb{C_n} \end{array}\end{align}

In the present work, the aerodynamic coefficient structure and the use of linear stability-derivative terms (e.g. dependence on $\alpha $ , $\beta $ and normalised rates $ {pb}/ {2{V_{air}}}, {qc}/ {2{V_{air}}}, {rb}/ {2{V_{air}}} $ ) follow standard flight-dynamics and control-oriented modeling practice [Reference Etkin and Reid26, Reference Stevens, Lewis and Johnson27]. Because a dedicated, high-fidelity aerodynamic database for the full F-35B STOVL envelope is not employed at this stage, the numerical values of the coefficient set ( ${C_{L0}}$ , ${C_{L{\rm{\alpha }}}}$ , ${C_{D0}}$ , ${K_{CD}}$ , ${C_{Y{\rm{\beta }}}}$ , ${C_{l{\rm{\beta }}}}$ , ${C_{lp}}$ , ${C_{lr}}$ , ${C_{m0}}$ , ${C_{m{\rm{\alpha }}}}$ , ${C_{mq}}$ , ${C_{n{\rm{\beta }}}}$ , ${C_{np}}$ , ${C_{nr}}$ ) were initialised from typical magnitudes used in control-oriented fixed-wing models and then refined to satisfy trim consistency at the selected approach condition and bounded/physically plausible force/moment levels relative to the vehicle weight and actuator limits. The aerodynamic coefficients used in this study should therefore be interpreted as representative control-oriented parameters rather than validated F-35B STOVL aerodynamic data. Replacing these coefficients with wind-tunnel, CFD, DATCOM or flight-test-identified values would improve the fidelity only within the validity range of the corresponding aerodynamic model structure. For the complete vertical/near-hover STOVL regime, however, a dedicated aerodynamic database including jet-induced flow, ground effect and multi-directional inflow would be required. The present study does not attempt to construct such a database; instead, it quantifies the magnitude of the resulting aerodynamic loads and uses them as bounded surrogate contributions in the integrated closed-loop simulation [Reference Williams and Vukelich28–Reference Blake30].

Since a full publicly available inertia tensor for the F-35B is not provided in open literature, the inertia parameters were initialised using control-oriented order-of-magnitude estimates and subsequently treated as tunable quantities.

The aircraft mass $m$ was selected as a nominal approach/landing value within the published weight class for the F-35B and may be adjusted to represent different payload/fuel conditions. Geometric parameters (wing area $S$ and wingspan $b$ ) were taken from publicly available F-35B aircraft specification sheets [31]. Mean aerodynamic chord $c$ enters the linear aerodynamic moment model mainly through the standard non-dimensional rate terms (e.g. ${qc}/{2{V_{air}}}$ ), which set the correct scaling of pitch-rate damping with respect to the aircraft’s reference length. When the detailed planform $c( y )$ is unavailable, we approximate the wing by an equivalent rectangular planform with the same reference area $S$ and span $b$ . For such a planform, the chord is $c = {S/b} $ providing a consistent order-of-magnitude reference length and yielding physically reasonable normalised-rate effects in a control-oriented simulation.

Aerodynamic coefficients were adopted as a low-speed linear surrogate model suitable for controller evaluation; rate-normalisation safeguards and angle saturations were included to prevent numerical singularities at very low airspeed. The parameters used in the aerodynamic force and moment model are listed in Table 2.

Table 2.

Representative control-oriented aerodynamic coefficients used in the low-order surrogate model

Table 2. Long description

A table with 12 rows and 4 columns. The columns are labeled Symbol/variable, Description, Recommended value, and Unit. The table lists various aerodynamic coefficients used in a low-order surrogate model. Row 1: Symbol/variable, Lift at zero angle-of-attack; Description, Lift at zero angle-of-attack; Recommended value, 0.30; Unit, blank. Row 2: Symbol/variable, Lift slope; Description, Lift slope; Recommended value, 4.5; Unit, 1/rad. Row 3: Symbol/variable, Zero-lift drag; Description, Zero-lift drag; Recommended value, 0.05; Unit, blank. Row 4: Symbol/variable, Induced drag factor; Description, Induced drag factor; Recommended value, 0.30; Unit, blank. Row 5: Symbol/variable, Side-force slope; Description, Side-force slope; Recommended value, 0.98; Unit, 1/rad. Row 6: Symbol/variable, Roll due to sideslip; Description, Roll due to sideslip; Recommended value, 0.12; Unit, 1/rad. Row 7: Symbol/variable, Roll damping; Description, Roll damping; Recommended value, 0.50; Unit, blank. Row 8: Symbol/variable, Roll due to yaw rate; Description, Roll due to yaw rate; Recommended value, 0.25; Unit, blank. Row 9: Symbol/variable, Pitch bias; Description, Pitch bias; Recommended value, 0.05; Unit, blank. Row 10: Symbol/variable, Pitch stability; Description, Pitch stability; Recommended value, 1.00; Unit, 1/rad. Row 11: Symbol/variable, Pitch damping; Description, Pitch damping; Recommended value, 8.0; Unit, blank. Row 12: Symbol/variable, Yaw due to sideslip; Description, Yaw due to sideslip; Recommended value, 0.25; Unit, 1/rad. Row 13: Symbol/variable, Yaw due to roll rate; Description, Yaw due to roll rate; Recommended value, 0.10; Unit, blank. Row 14: Symbol/variable, Yaw damping; Description, Yaw damping; Recommended value, 0.20; Unit, blank.

The aerodynamic model is retained in the present framework for a limited and explicitly control-oriented purpose. It is not intended to provide a validated prediction of the complete F-35B STOVL hover or near-ground aerodynamic environment. In the simulated low-speed inclined landing manoeuvre, the primary force and moment generation is provided by the lift fan, rear swivel nozzle and wing roll-post thrust channels, whereas no aerodynamic control-surface actuation is used by the controller. Therefore, the aerodynamic model acts only as a bounded passive airframe-load contribution in the nonlinear 6-DOF equations of motion. This avoids an entirely force-free aerodynamic assumption while keeping the model modular, so that the same simulation structure can later accept higher-fidelity aerodynamic data from wind-tunnel tests, CFD, DATCOM-type estimation or flight-test identification. In this sense, the aerodynamic model should be interpreted as a secondary surrogate load model included for completeness and extensibility, not as the dominant mechanism governing the landing response.

The aerodynamic load magnitudes reported here were obtained by post-processing the simulation outputs as ${F_{aero,mag}}(t) = \sqrt {F_{aero,x}^2 + F_{aero,y}^2 + F_{aero,z}^2} $ , ${M_{aero,mag}}(t) = \sqrt {M_{aero,x}^2 + M_{aero,y}^2 + M_{aero,z}^2} $ . The reported percentages were calculated using $100{F_{aero,max}}/\left( {mg} \right)$ , where ${F_{aero,max}} = {\mathrm{max}}\!\left( {{F_{aero,mag}}( t )} \right)$ . For the nominal landing case, a post-processing check showed that the maximum airspeed is 6.905 m/s, the maximum dynamic pressure is only 29.203 Pa, and the maximum aerodynamic force magnitude is 2,706.16 N, corresponding to approximately 1.379% of the aircraft weight. The maximum aerodynamic moment magnitude is 3,345.44 Nm. These values confirm that, in the present near-hover landing scenario, the aerodynamic loads remain secondary compared with the thrust and landing-gear contact loads. Therefore, the present aerodynamic formulation is used to provide bounded control-oriented force and moment contributions, while its validity remains limited at very low dynamic pressure, high angle-of-attack, jet-induced flow, ground-effect and multi-directional inflow conditions. The simulation results should thus be interpreted as a control-oriented integrated flight-dynamics and landing-gear assessment, not as a certified prediction of the full F-35B STOVL aerodynamic environment.

Landing behaviour is primarily driven by thrust-vectoring and landing-gear dynamics. The aerodynamic model does not claim high-fidelity prediction; it is retained merely as a transparent, bounded secondary surrogate for the 6-DOF simulation with explicitly acknowledged limitations.

2.4. Landing gear and suspension dynamics

To complete the physical representation of the F-35B, particularly for the critical transition from flight to ground contact, the rigid-body model is extended with a control-oriented landing gear and suspension subsystem. The aircraft is equipped with a simplified three-point landing gear model consisting of one front wheel and two rear wheels (left and right). This section details the interaction between the airframe and the landing surface, accounting for the mechanical properties of the struts and the kinematic constraints of the landing gear configuration. Figure 2 shows the schematic representation of the landing gear and suspension model.

Figure 2.

Schematic representation of the landing gear and suspension model. Base aircraft image adapted from Ref. (24); landing-gear spring–damper schematics and annotations were added by the authors.

The following assumptions are adopted for a control-oriented yet dynamically meaningful model:

• • The gear dynamics are dominated by vertical compliance (heave/roll/pitch coupling arises through attachment-point kinematics).

• • Suspension springs and dampers are linear in relative displacement/velocity.

• • Tire compliance is represented by a unilateral vertical spring; no tensile (pulling) force is allowed.

• • Runway is locally flat, but can be generalised via local runway heights.

• • Lateral tire forces (braking/side force) are not included in this section; only vertical support forces are modeled.

• • Longitudinal tire forces, wheel rotational dynamics, braking torque and rolling resistance are not included in the present model. The landing-gear subsystem is therefore used to compute vertical tire normal reactions, strut forces and suspension strokes, but it should not be interpreted as a complete ground-roll or wheel-spin model.

The landing gear system is modeled as a tricycle configuration consisting of one front wheel (index $f$ ) and two rear wheels (indices $w1$ for left and $w2$ for right). This subsystem introduces three additional degrees of freedom to the aircraft model, representing the vertical displacement of the unsprung masses ( ${z_{t,f}}$ , ${z_{t,w1}}$ , ${z_{t,w2}}$ ). Unlike the aircraft’s position in the NED frame ( $Z$ ), these coordinates are defined as upward-positive, following standard vehicle dynamics conventions. The circled single-leg schematics in Fig. 2 also illustrate the geometric meaning of the tire-compression variable. For the $i$ -th wheel, ${z_{t,i}}$ denotes the upward-positive unsprung wheel coordinate and ${z_{r,i}}$ denotes the local runway height under the same wheel. The tire compression is therefore defined as ${\delta _{t,i}} = {\textrm{max}}\!\left( {0,{z_{r,i}} - {z_{t,i}}} \right)$ . With this unilateral definition, ${\delta _{t,i}} = 0$ corresponds to an airborne wheel or no tire compression, and a tire-runway normal contact force is generated only when the wheel penetrates the local runway level in the simplified vertical contact model.

To determine the forces generated by the suspension, we must first calculate the instantaneous height and velocity of the attachment points (sprung points) on the airframe. The position of each gear attachment point in the inertial frame ( ${r_{g,iN}}$ ) is found by rotating the body-fixed offset vectors using the direction cosine matrix ( $\boldsymbol{R}_b^i$ ) defined in Section 2.1:

(71) \begin{align}\begin{array}{*{20}{c}}{\boldsymbol{r}_{g,iN}} = \boldsymbol{r} + \boldsymbol{R}_b^i{\boldsymbol{r}_{g,i}} \end{array}\end{align}

All displacement vectors are defined in Equations (72)–(74):

(72) \begin{align}\begin{array}{*{20}{c}}{{\boldsymbol{r}}_{g,f}} = \left[ {\begin{array}{*{20}{c}}{{L_g}}\\0\\{{L_c}}\end{array}} \right] \end{array}\end{align}

(73) \begin{align}\begin{array}{*{20}{c}}{{\boldsymbol{r}}_{g,w1}} = \left[ {\begin{array}{*{20}{c}}{ - {L_d}}\\{{L_e}}\\{{L_c}}\end{array}} \right] \end{array}\end{align}

(74) \begin{align}\begin{array}{*{20}{c}}{{\boldsymbol{r}}_{g,w2}} = \left[ {\begin{array}{*{20}{c}}{ - {L_d}}\\{{ - L_e}}\\{{L_c}}\end{array}} \right] \end{array}\end{align}

The vertical coordinate of the sprung point ( ${z_{s,i}}$ ) is derived by mapping the NED $Z$ components to the upward-positive suspension frame:

(75) \begin{align}\begin{array}{*{20}{c}}{z_{s,i}} = - {{\left( {{\boldsymbol{r}_{g,iN}}} \right)}_3} \end{array}\end{align}

Here ${(\!\cdot\!)}_3$ denotes the third component of a vector, i.e. the vertical (Down) component in the NED frame.

The velocity of the attachment point $i$ can be written using rigid-body point kinematics. Let ${{\textbf{v}}_b}$ be the CG velocity in the body frame and ${\boldsymbol{\omega} _b}$ the body angular rate. The point in the body frame is:

(76) \begin{align}\begin{array}{*{20}{c}}{{\textbf{v}}_{g,iB}} = {{\textbf{v}}_b} + {\boldsymbol{\omega} _b} \times {\boldsymbol{r}_{g,i}} \end{array}\end{align}

Using the DCM mapping of vectors:

\begin{align*}{{\textbf{v}}_{g,iN}} = \boldsymbol{R}_b^i{{\textbf{v}}_{g,iB}}\end{align*}

where ${{\textbf{v}}_{g,iN}}$ is the inertial velocity of the attachment point, accounting for the aircraft’s translational velocity and the angular velocity cross-product. The time derivative of the sprung point ( ${\dot z_{s,i}}$ ) is defined as:

(77) \begin{align}\begin{array}{*{20}{c}}{{\dot z}_{s,i}} = - {{\left( {{{\textbf{v}}_{g,iN}}} \right)}_3} \end{array}\end{align}

Each wheel assembly is treated as an independent dynamical system. The motion of the unsprung mass is governed by the balance between suspension force and tire reaction force. Using the unilateral tire-compression definition introduced above, the tire force for the $i$ -th wheel is written as:

(78) \begin{align}\begin{array}{*{20}{c}}{F_{tire,i}} = {k_{t,i\;}}{\delta _{t,i}} \end{array}\end{align}

Thus, the tire force is generated only when ${\delta _{t,i}} \gt 0$ ; otherwise, the wheel is airborne or unloaded and ${F_{tire,i}} = 0$ . During the contact phase, where ${z_{t,i}} \lt {z_{r,i}}$ Equation (78) reduces to:

(79) \begin{align}\begin{array}{*{20}{c}}\;{F_{tire,i}} = \;{k_{t,i\;}}\left( {{z_{r,i}}- {z_{t,i}}} \right) \end{array}\end{align}

Upon contact with the runway, the suspension force is modelled as a linear spring-damper system:

(80) \begin{align}\begin{array}{*{20}{c}}{F_{susp,i\;}} = \;{k_{s,i\;}}\left( {{z_{s,i}} - {z_{t,i}}} \right) + {c_{s,i\;}}\left( {{{\dot z}_{s,i}} - {{\dot z}_{t,i}}} \right)\; \end{array}\end{align}

For each wheel $i \in \left\{ {\,f,w1,w2} \right\}$ the vertical suspension force ( ${F_{susp,i\;}}$ ) and the vertical tire force ( ${F_{tire,i}}$ ) denote positive compression-force magnitudes. With the adopted sign convention for the unsprung-mass coordinate, the resulting dynamics is written as:

(81) \begin{equation} m_{t,i}\,\ddot{z}_{t,i} = -F_{susp,i} - F_{tire,i} \end{equation}

Substituting Equations (79) and (80) into (81) gives a closed-form of second-order system equation:

(82) \begin{equation} m_{t,i}\,\ddot{z}_{t,i} = - k_{s,i}\left(z_{s,i}-z_{t,i}\right) - c_{s,i}\left(\dot{z}_{s,i}-\dot{z}_{t,i}\right) - k_{t,i}\left(z_{r,i}-z_{t,i}\right) \end{equation}

This shows that wheel accelerations are driven by the relative motion between the airframe attachment point and the wheel, and by tire contact.

The forces generated by the suspension struts act as external loads on the rigid-body aircraft. Since these forces are calculated in the vertical suspension frame (upward-positive), they are first mapped into the NED frame as downward-positive forces:

(83) \begin{align}\begin{array}{*{20}{c}}{\boldsymbol{F}_{g,iN}} = \left[ {\begin{array}{*{20}{c}}0\\0\\{{F_{susp,i\;}}}\end{array}} \right] \end{array}\end{align}

These inertial forces are then rotated back into the aircraft’s body frame using the inverse rotation matrix $\boldsymbol{R}_n^i = {\left( {\boldsymbol{R}_b^i} \right)^{\mathrm{T}}}$ :

(84) \begin{align}\begin{array}{*{20}{c}}{\boldsymbol{F}_{g,i}} = \;\boldsymbol{R}_n^i{\boldsymbol{F}_{g,iN}} \end{array}\end{align}

The total gear force ( ${\boldsymbol{F}_{gear}}$ ) and the total moment ( ${\boldsymbol{M}_{gear}}$ ) about the centre of gravity are obtained by summing the contributions from all three wheels:

(85) \begin{align}\begin{array}{*{20}{c}}{\boldsymbol{F}_{gear}} = \;\mathop \sum \limits_{i \in \left\{ {\,f,w1,w2} \right\}} {\boldsymbol{F}_{g,i}} \end{array}\end{align}

(86) \begin{align}\begin{array}{*{20}{c}}{\boldsymbol{M}_{gear}} = \;\mathop \sum \limits_{i \in \left\{ {\,f,w1,w2} \right\}} {\boldsymbol{r}_{g,i}} \times {\boldsymbol{F}_{g,i}} \end{array}\end{align}

These vectors are then integrated into the governing equations of motion derived in Section 2.2, completing the 6-DOF dynamic loop for both flight and ground operations. Table 3 summarises the landing-gear and suspension parameters employed in the post-touchdown simulation.

Table 3.

Representative landing-gear and suspension parameters used in the control-oriented simulation model

Table 3. Long description

A table with 18 rows and 4 columns detailing landing-gear and suspension parameters used in a control-oriented simulation model. The columns are labeled Symbol/variable, Description, Representative value used, and Unit. The table includes various parameters such as front unsprung mass, rear-left unsprung mass, rear-right unsprung mass, front suspension stiffness, front suspension damping coefficient, rear-left suspension stiffness, rear-left suspension damping coefficient, rear-right suspension stiffness, rear-right suspension damping coefficient, front tire vertical stiffness, rear-left tire vertical stiffness, rear-right tire vertical stiffness, longitudinal offset from CG to front gear, longitudinal offset from CG to main gears, lateral half-track of main gears, vertical offset of each gear attachment point, assumed maximum allowable suspension stroke for each gear leg, and assumed maximum tire compression. Each row provides the representative value used and the corresponding unit for each parameter.

The landing-gear and suspension parameters listed in Table 3 are representative values used for the F-35B-type control-oriented simulation model. They should not be interpreted as certified, manufacturer-provided or flight-test-identified F-35B landing-gear design data. Since detailed F-35B strut, tire, damping and unsprung-mass properties are not publicly available, the parameters were selected as physically plausible order-of-magnitude values for a heavy STOVL aircraft landing simulation and were checked through the resulting touchdown loads and suspension-stroke responses. The equivalent uncompressed gear-leg reference length is taken as ${L_c} = 0.5$ m for the front, rear-left, and rear-right gears, consistent with the vertical attachment-point offset used in the simplified geometry. The allowable suspension compression stroke is assumed as ${\delta _{s,allow}} = 0.3$ m for each gear leg. These assumptions provide a transparent parameter basis for the present control-oriented simulation and can be replaced by validated aircraft-specific landing-gear data when such data become available.

5.1 Effect of landing-trajectory inclination angle

To further investigate the influence of the inclined landing geometry, an additional parametric study was conducted by varying the final horizontal landing distance ${X_f}$ while keeping the initial altitude ${h_0}$ , the reference manoeuvre duration ${T_m}$ , the total simulation time ${T_{sim}}$ , the controller gains, aircraft parameters and landing-gear parameters unchanged. The landing-trajectory inclination angle was defined as

(115) \begin{align}\begin{array}{*{20}{c}}\gamma = ta{n^{ - 1}}\left( {\dfrac{{{h_0}}}{{{X_f}}}} \right) \end{array}\end{align}

where ${h_0}$ is the initial altitude and ${X_f}$ is the prescribed final horizontal landing distance. This setup allows the effect of the landing-trajectory inclination to be evaluated separately from changes in manoeuvre duration, controller tuning or landing-gear parameters. The resulting cases are summarised in Table 5. The purpose of these cases is not to provide a full validation of the complete STOVL landing envelope, but to provide a controlled comparative assessment of the proposed closed-loop simulation framework under changes in the prescribed landing-trajectory geometry. In all three cases, the controller gain, aircraft parameters, landing-gear parameters, manoeuvre duration and simulation settings are kept unchanged; only the final horizontal landing distance ( ${X_f}$ ), and therefore the corresponding inclination angle, is varied.

Table 5.

Inclined landing cases used for the landing-trajectory inclination analysis

Table 5. Long description

A table comparing landing trajectory metrics for different inclined landing cases. The table has 11 rows and 4 columns. The columns are labeled Metric, Case A, Case B, and Case C. The row labels include Landing inclination angle, Maximum vertical tracking error, Maximum pitch angle, Maximum pitch rate, Maximum total thrust, Maximum nozzle angle, Peak tire normal load among all gears, Peak strut force among all gears, Max. suspension stroke among all gears, Stroke usage, Touchdown time of front gear, Touchdown time of rear-left gear, and Touchdown time of rear-right gear. Each case corresponds to a different final horizontal landing distance: 35 meters for Case A, 30 meters for Case B, and 25 meters for Case C. The values in the table represent various metrics measured during the landing trajectory analysis. For example, the Landing inclination angle is 23.199 degrees for Case A, 26.565 degrees for Case B, and 30.964 degrees for Case C. The Maximum vertical tracking error is 1.7334 meters for Case A, 1.5159 meters for Case B, and 1.3036 meters for Case C. The table provides a detailed comparison of these metrics across the three cases, allowing for an evaluation of the influence of landing-trajectory inclination.

The nominal case corresponds to ${X_f} = 30$ m and is the case used for the detailed time-history results in Figs 4–13. While Table 4 reports the individual front, rear-left and rear-right landing-gear responses for the nominal case, Table 6 uses the maximum value among the three landing gears as a compact metric for comparing different landing-trajectory inclination cases.

To visualise the geometric difference among the investigated landing-trajectory inclination cases, the reference and closed-loop landing trajectories are compared in Fig. 14 in the horizontal position–altitude plane. Here, the altitude is defined as ( $h = - Z$ ) because the inertial z-axis follows the NED convention. All three cases start from the same initial altitude ( ${h_0} = 15\,\text{m}$ ) and reach the touchdown plane ( $h = 0$ ), while the prescribed final horizontal distance ( ${X_f}$ ) is varied from 35 m to 25 m. Therefore, Case A corresponds to the mildest trajectory, Case B to the nominal trajectory and Case C to the steepest trajectory. The figure shows that reducing ( ${X_f}$ ) increases the landing-trajectory inclination angle, while the closed-loop trajectories remain close to their corresponding reference profiles.

Table 6.

Effect of landing-trajectory inclination on closed-loop tracking, control effort, touchdown load and landing-gear response

Table 6. Long description

A table comparing landing-gear touchdown and suspension-stroke characteristics for the front gear, rear-left gear, and rear-right gear. The table has 8 rows and 4 columns. The columns are labeled Quantity, Front gear, Rear-left gear, and Rear-right gear. The rows are labeled with different quantities and their respective values. Row 1: Touchdown time, [s], 7.99, 7.35, 7.35. Row 2: Maximum suspension stroke, [m], 0.0115, 0.0164, 0.0164. Row 3: Settled stroke after touchdown [m], 0.0095, 0.0130, 0.0130. Row 4: Assumed allowable stroke [m], 0.300, 0.300, 0.300. Row 5: Stroke usage [%], 3.85, 5.47, 5.47. Row 6: Peak tire normal load [N], 12,859, 14,958, 14,958. Row 7: Peak strut force [N], 11,093, 13,032, 13,032.

Figure 13.

Contact-dependent incremental suspension compression strokes during touchdown and post-touchdown settling: (top) front landing gear, (middle) rear-left landing gear and (bottom) rear-right landing gear. The vertical dashed lines indicate the corresponding touchdown instants.

Figure 14.

Comparison of the reference and closed-loop landing trajectories for the investigated inclination cases in the horizontal position–altitude plane.

Figure 14 visually confirms the expected effect of changing ( ${X_f}$ ) on the landing-trajectory geometry. Since all cases start from the same initial vertical position and use the same manoeuvre duration, reducing ( ${X_f}$ ) from 35 m to 25 m produces a progressively steeper landing trajectory in the ( $X$ )-( $Z$ ) plane. The closed-loop responses remain close to their corresponding reference trajectories for all three cases, indicating that the same controller configuration can track the investigated inclination cases without retuning. After this trajectory-level comparison, the corresponding quantitative performance metrics are summarised in Table 6.

Note: For each inclination case, ${F_{t,max}}$ , ${F_{s,max}}$ and ${\delta _{s,max}}$ represent the maximum values among the front, rear-left and rear-right landing gears. The nominal case ( ${X_f}$ =30 m) corresponds to the detailed landing-gear results reported in Table 4.

The comparison in Table 6 shows that the closed-loop responses remain bounded for all investigated landing-trajectory inclination angles. In the present setup, the inclination angle ( $\gamma $ ) is increased by reducing the prescribed horizontal landing distance ( ${X_f}$ ), while keeping ( ${h_0}$ ), ( ${T_m}$ ), ( ${T_{sim}}$ ), the controller gains and the landing-gear parameters unchanged. Therefore, the steeper cases require less horizontal translation within the same manoeuvre duration. As a result, the required pitch motion, pitch rate, total thrust and nozzle-vectoring demand decrease moderately as ( $\gamma $ ) increases. The touchdown loads and landing-gear responses remain bounded, with the maximum suspension stroke varying only slightly from 0.016407 m to 0.015986 m over the investigated range. In all cases, the maximum stroke usage remains below 6% of the assumed allowable suspension stroke, indicating that the landing-gear subsystem operates within the prescribed working range.

To make the differences among the inclination cases more directly observable, the numerical comparison in Table 6 is complemented by the time-history responses shown in Fig. 15. The vertical tracking error and total thrust-command magnitude are selected as representative quantities because they describe the main closed-loop tracking performance and actuator-demand trends without duplicating the detailed nominal-case responses already presented in the preceding figures.

Figure 15.

Time-history comparison of the investigated landing-trajectory inclination cases: (top) vertical tracking error and (bottom) total thrust-command magnitude.

Figure 15 confirms the trends observed in Table 6. Case A produces the largest vertical tracking error and total thrust-command demand, whereas Case C gives the lowest values among the investigated cases. This trend is mainly associated with the larger longitudinal displacement required in Case A over the same manoeuvre duration. The responses remain bounded for all three cases, indicating that the same controller configuration can handle the investigated changes in landing-trajectory inclination without retuning. The landing-gear-related quantities are reported in Table 6, where the maximum suspension stroke remains well below the assumed allowable stroke for all cases. The inclination-angle cases are not intended to constitute a full validation of the complete STOVL landing envelope. Instead, they provide a controlled comparative assessment of the proposed closed-loop simulation framework under changes in the prescribed landing-trajectory geometry. In all three cases, the controller gain, aircraft parameters, landing-gear parameters, manoeuvre duration and simulation settings are kept unchanged. Therefore, the comparison is used to examine the resulting trends in tracking performance, actuator demand and landing-gear response under otherwise identical modelling assumptions.

To provide a compact visual interpretation of the inclination-angle study, the main performance indicators from Table 6 are plotted in Fig. 16 as functions of the landing-trajectory inclination angle ( $\gamma $ ). The selected metrics represent the four aspects emphasised in the parametric comparison: vertical tracking performance, control effort, touchdown loading and landing-gear compression. Specifically, Fig. 16 compares the maximum vertical tracking error, maximum total thrust, peak tire normal load among all gears and maximum suspension stroke among all gears for the three investigated inclination cases. Note that the markers indicate the simulated inclination cases, and the connecting lines are used only to visualise the trend.

Figure 16.

Effect of landing-trajectory inclination angle on representative performance metrics: (a) maximum vertical tracking error, (b) maximum total thrust, (c) peak tire normal load among all gears and (d) maximum suspension stroke among all gears.

The comparison in Fig. 16 supports the quantitative trends reported in Table 6. The maximum vertical tracking error decreases from 1.7334 m to 1.3036 m as the landing-trajectory inclination angle increases from ( ${23.199^\circ }$ ) to ( ${30.964^\circ }$ ). A similar decreasing trend is observed in the pitch-related control demand and actuator effort, as the maximum total thrust decreases from ( $2.4722 \times {10^5}$ ) N to ( $2.2132 \times {10^5}$ ) N and the maximum nozzle angle decreases from ( ${3.4210^\circ }$ ) to ( ${2.8712^\circ }$ ). This behaviour is consistent with the present parameterisation of the inclination study: increasing ( $\gamma $ ) is achieved by reducing ( ${X_f}$ ) while keeping ( ${h_0}$ ) and ( ${T_m}$ ) constant, so the steeper cases require less horizontal translation within the same manoeuvre duration. Consequently, the longitudinal pitch and nozzle-vectoring demand are reduced. The touchdown-related quantities also remain bounded over the investigated range. The peak tire normal load decreases from 15,539 N to 14,351 N, while the maximum suspension stroke changes only slightly from 0.016407 m to 0.015986 m. These results indicate that the proposed closed-loop model maintains bounded responses for all investigated inclination cases and that the landing-gear subsystem remains within the assumed working range.

To make the low-speed nature of the nominal inclined landing manoeuvre explicit, the velocity responses for the nominal case are shown in Fig. 17.

Figure 17.

Velocity responses for the nominal inclined landing case: (top) horizontal inertial velocity $\,\dot{\!X}$ , (middle) vertical inertial velocity $\,\dot{\!Z}$ under the NED convention and (bottom) total airspeed ${V_{air}}$ .

Figure 17 confirms the low-speed near-hover character of the nominal inclined landing manoeuvre. The horizontal inertial velocity $\,\dot{\!X}$ reaches a maximum value of approximately 6.0 m/s, which is consistent with the fifth-order reference trajectory for ${X_f} = 30$ m and ${T_m} = 10$ s, whose reference peak is approximately 5.63 m/s. Similarly, the vertical inertial velocity $\,\dot{\!Z}$ remains within a low-speed descent range, with a peak value of approximately 3.5 m/s, compared with the reference peak of approximately 2.81 m/s. The short negative excursion of $\,\dot{\!Z}$ near touchdown represents a transient upward correction under the NED sign convention during the contact and post-touchdown settling phase. The total airspeed reaches a maximum value of approximately 6.9 m/s, which is consistent with the post-processing value ${V_{air,max}} = 6.905$ m/s reported in Section 2.3. Therefore, the velocity responses support the interpretation that the simulated landing is a low-speed inclined touchdown manoeuvre and that the aerodynamic loads remain secondary compared with the thrust and landing-gear contact loads.

5.2 Wind-gust-equivalent external-force disturbance case

In this section, the disturbance case is treated as a gust-equivalent external-force perturbation rather than as a full air-relative wind-field simulation. Therefore, the disturbance is not introduced through the aerodynamic relative-velocity formulation of Section 2.3. Instead, bounded body-frame disturbance forces are directly applied to the nonlinear 6-DOF equations of motion in order to assess the closed-loop response of the controller and landing-gear subsystem under an additional transient load. This approach is useful for implementation-level disturbance assessment, but it should not be interpreted as a complete aerodynamic wind-gust model. A full wind-field simulation, especially for wind speeds comparable to the aircraft speed, would require a validated low-speed/STOVL aerodynamic model capable of representing multi-directional inflow, jet-induced flow and ground-effect interactions.

To provide an additional assessment beyond the nominal landing condition, a bounded wind-gust-equivalent disturbance case was considered for the nominal inclination case m ${X_f} = 30$ . Since the present study is primarily focused on a control-oriented integrated flight-dynamics and landing-gear simulation framework, the gust effect was introduced as an equivalent external disturbance force acting in the body frame through the existing disturbance-force channels of the translational equations. The disturbance was applied between $t = 3$ s and $t = 9$ s as

(116) \begin{align}\begin{array}{*{20}{c}}{F_{dist,x}}(t) = 0.03W\sin \left( {\dfrac{{2\pi \left( {t - 3} \right)}}{3}} \right) \end{array}\end{align}

(117) \begin{align}\begin{array}{*{20}{c}}{F_{dist,z}}(t) = 0.04W\sin \left( {\dfrac{{2\pi \left( {t - 3} \right)}}{2}} \right) \end{array}\end{align}

where $W = mg$ is the aircraft weight. The lateral disturbance force was set to zero in this preliminary disturbance case. The nominal case corresponds to ${F_{dist,x}} = {F_{dist,y}} = {F_{dist,z}} = 0$ . This simplified deterministic disturbance input is not intended to represent a full stochastic turbulence or ship-airwake model, but it provides a compact robustness-oriented check of the closed-loop response under bounded external forcing.

The selected disturbance amplitudes are used as an equivalent force-level representation of a transient gust/load perturbation. Since these forces are applied directly to the body-frame equations of motion, they do not modify the aerodynamic angles, dynamic pressure or air-relative velocity used in the low-order aerodynamic surrogate model. Consequently, the results in this section should be interpreted as a closed-loop response to an external force disturbance, not as a high-fidelity prediction of aircraft behaviour in a real atmospheric wind field.

Table 7 shows that the gust-equivalent disturbance mainly affects the attitude-rate response, actuator demand and landing-gear loads, while the maximum vertical tracking error remains nearly unchanged compared with the nominal case. The maximum total thrust increases from ( $2.2907 \times {10^5}$ ) N to ( $2.4326 \times {10^5}$ ) N, and the maximum pitch rate increases from 0.2082 rad/s to 0.26014 rad/s, corresponding to an increase of approximately 25%. The touchdown-related quantities are more sensitive to the disturbance: the peak tire normal load increases from 14958 N to 16086 N, the peak strut force increases from 13032 N to 32999 N and the maximum contact-dependent suspension stroke increases from 0.016380 m to 0.053848 m. Although the suspension stroke increases noticeably under the gust-equivalent disturbance, the corresponding stroke usage remains 17.95%, which is still well below the assumed allowable suspension stroke. Therefore, the additional disturbance case indicates that the proposed simulation framework can be used not only for nominal touchdown analysis, but also for preliminary disturbance-sensitivity assessment. A more complete treatment of stochastic turbulence, ship-airwake effects and robust or gain-scheduled controller synthesis is left for future work.

Table 7.

Effect of the wind-gust-equivalent disturbance on the nominal landing case and summary of the main response changes

Table 7. Long description

A table with seven rows and five columns comparing the effects of wind-gust-equivalent disturbance on various metrics of a nominal landing case. The columns are labeled ‘Metric’, ‘Nominal case (Case B)’, ‘Gust-equivalent disturbance case’, and ‘Main observation’. The rows list different metrics and their values under both nominal and gust-equivalent disturbance cases, along with observations of the changes. The metrics include maximum vertical tracking error, maximum pitch angle, maximum pitch rate, maximum total thrust, maximum nozzle angle, peak tire normal load among all gears, peak strut force among all gears, maximum suspension stroke among all gears, and stroke usage. Each row provides specific values for these metrics under both conditions and notes the main observations, such as slight increases, significant increases, or nearly unchanged values. Notable trends include significant increases in peak strut force and maximum suspension stroke under the gust-equivalent disturbance case, while the maximum vertical tracking error remains nearly unchanged.

To make the effect of the gust-equivalent external-force disturbance more directly observable, Fig. 18 compares representative time histories of the nominal and disturbed cases. The pitch angle is selected to illustrate the aircraft attitude response, while the maximum contact-dependent suspension stroke is selected to show the post-touchdown landing-gear response, where the disturbance effect is more pronounced.

Figure 18.

Time-history comparison of the nominal case and the gust-equivalent external-force disturbance case: (top) pitch angle $\theta $ , (bottom) maximum contact-dependent suspension stroke ${\delta _{s,max}}$ .

Figure 18 shows that the gust-equivalent external-force disturbance causes only a modest change in the pitch-angle response but produces a more pronounced increase in the post-touchdown suspension compression response. The maximum contact-dependent suspension stroke increases from 0.0164 m in the nominal case to 0.0538 m in the disturbed case. Although this corresponds to a noticeable increase in the landing-gear response, the maximum value remains well below the assumed allowable suspension stroke of 0.3 m. Therefore, the disturbance mainly increases the touchdown/settling loads and suspension motion, while the overall closed-loop response remains bounded. Since the disturbance is applied as an external body-frame force rather than through a full air-relative wind-field model, these results should be interpreted as an implementation-level disturbance assessment rather than a complete aerodynamic wind-gust validation.