Free flight simulation in OpenVOGEL

One of the most recent features developed in the OpenVOGEL suit is the free flight simulation module. This module combines the unsteady aeroelastic solver (which provides air loads for a given flow field) with the numerical integration of the rigid body equations of motion (in 6 degrees of freedom). It resembles in some aspects the aeroelastic module, but here the flow field is directly manipulated from the outcome of the equations (in the aeroelastic module, the variation of the flow field is a consequence of the repositioning if the shedding edges).

Coupling algorithm

The free flight simulation consists of an algorithm that couples the air loads and the motion. The air-loads are first computed using the instantaneous flow field at the beginning of each time step using the previous state. Based on these air loads, the motion is predicted (through numerical integration) for the next time step and the result is translated into a new flow field. The air loads are then recomputed in the new state and the motion is corrected in a cycle that is repeated until convergence is reached. Once the convergence criteria are met, the time is advanced by one step.

Equations of motion

If we attempt to write the equations of motion of a rigid body in an inertial reference frame, the inertial properties of the body become time dependent. So to reduce the complexity of the equations, the motion can better be described from the point of view of the body itself (a non inertial reference frame), where the inertia tensor is an invariant. Additionally, to avoid having the kinetic moment being coupled to the momentum, it is convenient to take the reference point for the moments as the center of gravity. Finally, to avoid a complex coupling of the angular velocities in the angular momentum equations, it is much simpler to align the axes with the main axes of inertia.

With all this said, the equations of motion for the rigid aircraft model including the external forces and gravity field can be written as follows ^[1]:

${\frac {^{R}d}{dt}}{\mathbf {G}}={\mathbf {F}}+m{\mathbf {g}}$

${\frac {^{R}d}{dt}}{\mathbf {H}}={\mathbf {M}}$

Note that we have used Tenembaum's notation, in which the superscript R at the left of the differential operator means that the derivatives are taken with respect to the inertial reference frame (for the sake of simplicity, the earth in our case).

To convert the derivatives to the non inertial reference frame A (the moving aircraft), we use the reference frame transformation formula:

${\frac {^{R}d}{dt}}{\mathbf {a}}={\frac {^{A}d}{dt}}{\mathbf {a}}+{\mathbf {\omega }}\times {\mathbf {a}}$

where ${\mathbf {a}}$ is a generic vector.

The equations turn now to be:

${\frac {^{A}d}{dt}}{\mathbf {v}}+{\mathbf {\omega }}\times {\mathbf {v}}={\frac {1}{m}}{\mathbf {F}}+{\mathbf {g}}$

$II.{\frac {^{A}d}{dt}}{\mathbf {\omega }}+{\mathbf {\omega }}\times II.{\mathbf {\omega }}={\mathbf {M}}$

Tenembaum would accurately write the angular velocity vector as ${^{R}}{\mathbf {\omega }}^{A}$ , however, to simplify the notation (and because we only have two reference frames here) we simply write ${\mathbf {\omega }}$ to refer to the angular velocity of the aircraft with respect to the earth. The matrix $II$ is the tensor of inertia, which for us is an identity matrix with $I_{xx}$ , $I_{yy}$ and $I_{zz}$ in the main diagonal and 0's anywhere else.

In order to use the equations for numerical computations, it is necessary to write them in scalar form and to leave the derivatives alone at the left side:

Momentum

 ${\frac {^{A}d}{dt}}v_{x}={\frac {F_{x}}{m}}+g_{x}-\omega _{y}v_{z}+\omega _{z}v_{y}$ 

 ${\frac {^{A}d}{dt}}v_{y}={\frac {F_{y}}{m}}+g_{y}-\omega _{z}v_{x}+\omega _{x}v_{z}$ 

 ${\frac {^{A}d}{dt}}v_{z}={\frac {F_{z}}{m}}+g_{z}-\omega _{x}v_{y}+\omega _{y}v_{x}$ 

Angular momentum (Euler's equations)

 ${\frac {^{A}d}{dt}}\omega _{x}={\frac {1}{I_{xx}}}[M_{x}+(I_{yy}-I_{zz})\omega _{y}\omega _{z}]$ 

 ${\frac {^{A}d}{dt}}\omega _{y}={\frac {1}{I_{yy}}}[M_{y}+(I_{zz}-I_{xx})\omega _{z}\omega _{x}]$ 

 ${\frac {^{A}d}{dt}}\omega _{z}={\frac {1}{I_{zz}}}[M_{z}+(I_{xx}-I_{yy})\omega _{x}\omega _{y}]$

These six equations are not sufficient to describe the motion, since the gravity acceleration vector ${\mathbf {g}}$ also varies with time in the non inertial reference frame A. So we need to add the next three kinematic equations to the set:

${\frac {^{R}d}{dt}}{\mathbf {g}}={\mathbf {0}}\Rightarrow {\frac {^{A}d}{dt}}{\mathbf {g}}=-{\mathbf {\omega }}\times {\mathbf {g}}$

and thus, in scalar form:

Rotation of the gravity field

 ${\frac {^{A}d}{dt}}g_{x}=\omega _{z}g_{y}-\omega _{y}g_{z}$ 

 ${\frac {^{A}d}{dt}}g_{y}=\omega _{x}g_{z}-\omega _{z}g_{x}$ 

 ${\frac {^{A}d}{dt}}g_{z}=\omega _{y}g_{x}-\omega _{x}g_{y}$

Finally, if we want to output the position of the aircraft, then we need to integrate it in parallel also:

Position

 ${\frac {^{R}d}{dt}}P_{x}=v_{x}+\omega _{y}v_{z}-\omega _{z}v_{y}$ 

 ${\frac {^{R}d}{dt}}P_{y}=v_{y}+\omega _{z}v_{x}-\omega _{x}v_{z}$ 

 ${\frac {^{R}d}{dt}}P_{z}=v_{z}+\omega _{x}v_{y}-\omega _{y}v_{x}$

Note that this position will be relative to the aircraft reference coordinates, so to make it mean full, all points should be remapped to the global reference frame. The program does this when loading the results using the last orientation of the model.

Numerical integration

All the equations derived above form an initial value problem with 12 variables that can be solved using an appropriate numerical integration algorithm. One would probably think of the Runge-Kutta methods as first option, however they are not a good option since they require evaluating the forces at an intermediate time step, creating a complex scheme in our case (the UVLM method is quite heavy to be called continuously). Therefore, probably the best option is to use some kind of predictor corrector method that uses a fixed time step.

The numerical algorithm programmed in OpenVOGEL follows an scheme from Preidikman. This scheme is self-starting. It starts with the Euler method, then it switches to two Adams-Bashford and Adams-Moulton steps of increasing grade, and from the fourth step on it uses the Hamming method only.

In the next frames, vector X is the array (or better said, the record) containing all independent variables of the equations, and vector DX is the array containing the associated derivatives.

Time step 1

> Predict using initial derivatives (Euler)

 X(1) = X(0) + Dt * DX(0)

> Correct for K=0 to M (Modified Euler)

 Calculate new forces and derivatives...

 X(1) = X(0) + (0.5# * Dt) * (DX(0) + DX(1))

Time step 2

> Predict using previous derivatives (Adams-Bashford)

 X(2) = X(1) + (0.5# * Dt) * (3.0# * DX(1) - DX(0))

> Correct for K=0 to M (Adams-Moulton)

 Calculate new forces and derivatives...

 X(2) = X(1) + (Dt / 12.0#) * (5.0# * DX(2) + 8.0# * DX(1) - DX(0))

Time step 3

> Predict using previous derivatives (Adams-Bashford)

 X(3) = X(2) + (Dt / 12.0#) * (23.0# * DX(2) - 16.0# * DX(1) + 5.0# * DX(0))

> Correct for K=0 to M (Adams-Moulton)

 Calculate new forces and derivatives...

 X(3) = X(2) + Dt / 24.0# * (9.0# * DX(3) + 19.0# * DX(2) - 5.0# * DX(1) + DX(0))

 TE(3) = (9.0# / 121.0#) * (X(3) - XS)

Time steps S = 4 to N (Hamming)

> Predict using previous derivatives

 XP = X(S - 4) + (Dt * 4.0# / 3.0#) * (2.0 * DX(S - 1) - DX(S - 2) + 2.0 * DX(S - 3))
 
 X(S) = XP + 112.0# / 9.0# * TE(S - 1)

> Correct for K=0 to M

 Calculate new forces and derivatives...

 X(S) = (1.0# / 8.0#) * (9.0# * X(S - 1) - X(S - 3) + 3.0# * Dt * (DX(S) + 2.0# * DX(S - 1) - DX(S - 2)))
 
 TE(S) = (9.0# / 121.0#) * ((X(S) - XP))
 
 X(S) = X(S) - TE(S)

Validation

Our VB.NET implementation of the algorithm has been validated for a simple 1-D harmonic oscillator. In such a system the force is a linear function of the position and the velocity:

$F_{x}=-K.P_{x}-C.v_{x}$

The analytic solution for the displacement considering only an initial velocity is ^[2]^[3]:

$P_{x}(t)=e^{-\xi \omega t}v_{x}(0)sin(\omega _{d}t)$

where:

$\omega ={\sqrt {\frac {K}{M}}}$

$C_{c}=2{\sqrt {KM}}$

$\xi ={\frac {C}{C_{c}}}$

$\omega _{d}=\omega {\sqrt {1-\xi ^{2}}}$

The agreement of the numerical solution for this problem is excellent, as can be seen on the previous graphic. With 400 time steps and a time inteval of 0.025s, the simulation time was 10 seconds. Along the complete integration domain, the error could be maintained under 0.5%.

Reference

↑ "Fundamentals of applied dynamics", Roberto A. Tenembaum
↑ "Introducción a la dinámica estructural", Julio Massa
↑ "Física I", Marcelo Alonso, Edward J. Finn

[1] "Fundamentals of applied dynamics", Roberto A. Tenembaum

[2] "Introducción a la dinámica estructural", Julio Massa

[3] "Física I", Marcelo Alonso, Edward J. Finn

[1]

[2]

[3]