OpenVOGEL/自由飞行

OpenVOGEL 中的自由飞行模拟

OpenVOGEL 套件中最近开发的最重要的功能之一是自由飞行模拟模块。该模块将非定常气动弹性求解器（为给定流场提供气动力）与刚体运动方程（在 6 个自由度内）的数值积分相结合。在某些方面，它类似于气动弹性模块，但在这里，流场直接由方程的结果操纵（在气动弹性模块中，流场的变化是脱落边重新定位的结果）。

耦合算法

自由飞行模拟包括一个算法，该算法将气动力和运动耦合在一起。首先使用前一个状态在每个时间步长的开始使用瞬时流场计算气动力。基于这些气动力，预测下一时间步的运动（通过数值积分），并将结果转换为新的流场。然后在新的状态下重新计算气动力，并在循环中对运动进行校正，该循环会一直重复，直到达到收敛。一旦满足收敛条件，时间就会前进一步。

运动方程

如果我们尝试在惯性参考系中写出刚体的运动方程，那么物体的惯性属性就会随时间变化。因此，为了降低方程的复杂性，最好从物体本身的观点（非惯性参考系）来描述运动，在该参考系中，惯性张量是不变的。此外，为了避免动量矩与动量耦合，方便地将力矩的参考点取为重心。最后，为了避免角动量方程中角速度的复杂耦合，将轴与惯性的主要轴对齐要简单得多。

综上所述，包括外力和重力场在内的刚性飞机模型的运动方程可以写成如下形式

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

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

请注意，我们使用了Tenembaum 的符号，其中微分运算符左侧的R 上标表示导数是相对于惯性参考系（为了简单起见，我们这里采用地球）而言的。

为了将导数转换为非惯性参考系A（运动中的飞机），我们使用参考系转换公式

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

其中 ${\mathbf {a}}$ 是一个泛型向量。

现在方程变为

${\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 会将角速度矢量准确地写成 ${^{R}}{\mathbf {\omega }}^{A}$，但是为了简化符号（因为我们这里只有两个参考系），我们简单地写成 ${\mathbf {\omega }}$ 来表示飞机相对于地球的角速度。矩阵 $II$ 是惯性张量，对我们来说是一个单位矩阵，其主对角线上有 $I_{xx}$，$I_{yy}$ 和 $I_{zz}$，其他地方都是 0。

为了将方程用于数值计算，需要将它们写成标量形式，并将导数保留在左侧

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}]$

这六个方程不足以描述运动，因为重力加速度矢量 ${\mathbf {g}}$ 在非惯性参考系 A 中也随时间变化。因此，我们需要将接下来的三个运动学方程添加到方程组中

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

因此，以标量形式

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}$

最后，如果我们要输出飞机的位置，那么我们还需要并行地对其进行积分

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}$

请注意，此位置将相对于飞机参考坐标，因此为了使其具有完全意义，所有点都应重新映射到全局参考系。程序在使用模型的最后方位加载结果时执行此操作。

数值积分

上面推导的所有方程形成了一个具有 12 个变量的初值问题，可以使用适当的数值积分算法求解。人们可能会首先想到龙格-库塔方法，但它们不是一个好的选择，因为它们需要在中间时间步长评估力，在我们这种情况下会创建一个复杂的方案（UVLM 方法非常繁重，无法连续调用）。因此，最好的选择可能是使用某种使用固定时间步长的预估校正方法。

OpenVOGEL 中编写的数值算法遵循 Preidikman 的方案。该方案是自启动的。它从欧拉方法开始，然后切换到两个越来越高阶的阿达姆斯-巴什福斯和阿达姆斯-莫尔顿步骤，从第四步起，它仅使用汉明方法。

在接下来的帧中，向量 X 是包含方程所有自变量的数组（或者更确切地说，是记录），向量 DX 是包含相关导数的数组。

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)

验证

使用简单谐振子验证数值积分算法

我们对该算法的 VB.NET 实现已针对简单的 1-D 谐振子进行了验证。在这样的系统中，力是位置和速度的线性函数

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

仅考虑初始速度的位移的解析解为

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

其中

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

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

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

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

如先前图形所示，该问题的数值解吻合得非常好。使用 400 个时间步长和 0.025 秒的时间间隔，模拟时间为 10 秒。在整个积分域内，误差可以保持在 0.5% 以下。