## Subsystem Aircraft Equations of Motion (Beaver)

The dynamics of a rigid aircraft can be expressed by twelve Ordinary Differential Equations (ODEs). The equations which determine the magnitude of the time-derivatives of the twelve state variables as a function of the external forces and moments, external disturbances, and the states themselves, have been gathered in the subsystem Aircraft Equations of Motion (Beaver) of the Beaver model. This subsystem also contains the integrator which obtains the states from their time-derivatives.

The twelve ODEs themselves are independent of the aircraft considered. However, some of the aerodynamic forces or moments may depend upon the time-derivatives of the state variables, which would make the differential equations implicit. It depends upon the resulting implicit expressions whether or not Simulink can cope with them. However, since it is usually not very difficult to rewrite the implicit differential equation(s) in terms of a set generalized explicit ODEs and a separate aircraft-dependent correction term, an aircraft-dependent block xdotcorr was included in the subsystem Aircraft Equations of Motion, which explains why the full subsystem name has become Aircraft Equations of Motion (Beaver).

Aircraft Equations of Motion (Beaver) contains the following blocks:

- xdotcorr (Beaver) contains the aircraft-dependent corrections to the state equations which were needed to convert the original implicit differential equations into a set of explicit ODEs (obviously, this block is currently valid for the DHC-2 'Beaver' aircraft).
- uvw contains the transformation equations which yield the components of the true airspeed along the body-axes of the aircraft, given the true airspeed, angle of attack, and sideslip angle.
- an Integrator block is used to obtain the time-trajectories of the state variables themselves from their derivatives with respect to time.
- the Gain block XFIX makes it possible to artificially fix states, by performing an element-by-element multiplication of the vector xdot with the vector xfix. All elements of xfix either equal one or zero, meaning that the corresponding states are either free to vary or artificially kept constant; xfix can be set manually or defined by means of the utility FIXSTATE.
- 12 ODEs contains the twelve ODEs themselves.

### Computation order

First, the velocity components u, v, and w are computed, using the current value of the state vector x, and the outputvector yhlp from the masked subsystem Hlpfcn (this vector contains frequently used sines and cosines of angular state variables). The velocity components u, v, and w are needed for solving the differential equations concerning the aircraft-coordinates xe and ye. In order to account for contributions of non-steady atmosphere, the velocity components due to wind and/or atmospheric turbulence, uw, vw, and ww, are extracted from the vector uwind, and added to u, v, and w. Next, the time-derivatives of the states are computed in the subsystem 12 ODEs, which uses the current state vector x, the helpvector yhlp, the vectors Ftot and Mtot (containing the total forces and moments), and the vector [u+uw, v+vw, w+ww]T as inputs.

The time-derivatives of the states are gathered in the vector xdot = dx/dt. As said before, the vector xdot that leaves the subsystem 12 ODEs needs to be corrected for aircraft-dependent terms. This is done in the subsystem xdotcorr (Beaver), which needs the current state vector x, and, in the case of the 'Beaver' aircraft, the outputvector from the masked subsystem Atmosph in the subsystem Airdata Group (to be more specific: xdotcorr needs to know the air-density for the correction of the beta-dot equation). The outputvector from xdotcorr, the corrected version of xdot, is send to the gain XFIX, which is used to artificially fix states if required by the user. Finally, the time-derivatives are sent to the Integrator block, where they are integrated to obtain the updated value of the state vector x. Note: the Integrator block uses the vector xinco as initial condition (see below).

The state vector x is fed back to the other subsystems in Beaver, and also internally in the subsystem Aircraft Equations of Motion, because the time-derivatives of the state variables depend upon the state variables themselves. The second level of Beaver contains a masked subsystem block Hlpfcn, which computes the sines and cosines of the angles alpha, beta, psi, theta, and phi. These 'helpvariables' are also fed back to other subsystems in Beaver, including Aircraft Equations of Motion. The resulting feedback loops can be distinguished from the other signals by their magenta foreground color (see also the color definitions).

### Inputs to the subsystem Aircraft Equations of Motion (Beaver)

      Ftot  = [Fx Fy Fz]'   (total forces along body-axes, from FMsort)
      Mtot  = [L M N]'     (total moments along body-axes, from FMsort)
      uwind = [uw vw ww uwdot vwdot wwdot]'  (atmospheric disturbances)
      yatm  = [rho ps T mu g]'        (atmosphere variables, computed
                                       in the masked subsystem Atmosph)
      yhlp  = [cos(alpha) sin(alpha) cos(beta) sin(beta) tan(beta) ...
               ... sin(psi) cos(psi) sin(theta) cos(theta) sin(phi) ...
               ... cos(phi) ]'           (sines & cosines, from Hlpfcn)

      Fx, Fy, Fz: total external forces along body-axes [N]
      L, M, N   : total external moments along body-axes [Nm]
      uw, vw, ww: body-axes components of velocity component due to
                  non-steady atmosphere (wind + turbulence) [m/s]
     {uwdot     : d(uw)/dt [m/s^2]                                   }
     {vwdot     : d(vw)/dt [m/s^2]                                   }
     {wwdot     : d(ww)/dt [m/s^2]                                   }
      rho       : air density [kg/m^3]
     {ps        : static pressure [N/m^2]                            }
     {T         : air temperature [K]                                }
     {mu        : dynamic viscosity [kg/(m*s)]                       }
     {g         : acceleration of gravity [m/s^2]                    }

The variables which are not used by any of the (masked) subsystem blocks from the subsystem Aircraft Equations of Motion have been put between curly braces. In addition to these variables, all elements from yhlp are used.

### Outputs from the subsystem Aircraft Equations of Motion (Beaver)

      x     = [V alpha beta p q r psi theta phi xe ye H]'    (states)
      xdot  = dx/dt                      (time-derivatives of states)
      ybvel = [u v w]'        (body-axes components of true airspeed)

      V    : true airspeed [m/s]
      alpha: angle of attack [rad]
      beta : sideslip angle [rad]
      p    : roll-rate [rad/s]
      q    : pitch-rate [rad/s]
      r    : yaw-rate [rad/s]
      psi  : yaw-angle [rad]
      theta: pitch-angle [rad]
      phi  : roll-angle [rad]
      xe   : X-coordinate, relative to Earth-axes [m]
      ye   : Y-coordinate, relative to Earth-axes [m]
      H    : altitude above sea level [m]
      u    : component of V along XB-axis [m/s]
      v    : component of V along YB-axis [m/s]
      w    : component of V along ZB-axis [m/s]

### Parameters to be defined in the Matlab workspace

- xinco: Initial value of the state vector x. Usually, a steady-state trimmed-flight condition is used as initial condition for the aircraft ('condition' = combination of input variables and state vector, so xinco does not specify a complete initial flight condition by itself!). Trimmed-flight conditions can be computed with the trim-routine ACTRIM. Initial conditions can also be retrieved from file by means of the Matlab utility TRILOAD.
- The block Integrator specifies xinco as initial value of the state vector.
- xfix: Vector with twelve elements (length(xfix) = length(x)) which are all equal to either one or zero. If the nth element of xfix is equal to one, the corresponding element of x (= the corresponding state variable) may vary freely, but if the nth element of xfix is equal to zero, the corresponding state is artificially fixed to its initial value, no matter what its actual value should be. This may be useful for the analysis of some types of aircraft responses. The vector xfix can be defined manually, or by running FIXSTATE before starting a new simulation. Usually, xfix = ones(1,12) will be used for simulations. This default value will be set automatically when a new simulation is started while xfix has not yet been defined in the workspace.
- GM1: Matrix which contains some basic geometric properties of the aircraft, and its mass. GM1 is needed by the subsystem 12 ODEs.
- GM2: Matrix with inertia parameters of the aircraft. GM2 is needed by the subsystem 12 ODEs.