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% 以下。