1. Область модели
Самолет считается жестким телом. Аэродинамика считается по набору панелей:
крыло, хвостовое оперение и киль. Каждая панель видит собственный локальный поток,
получает свои \(C_L\), \(C_D\), \(C_m\), после чего силы и моменты суммируются на центре масс.
Это приближенная инженерная модель для малых планеров. Она не заменяет CFD,
X-Plane blade element model или JSBSim FDM, но дает прозрачные формулы для отладки.
2. Константы и единицы
| Величина | Значение | Комментарий |
|---|
| \(\rho\) | \(1.225\ kg/m^3\) | плотность воздуха у земли |
| \(g\) | \(9.81\ m/s^2\) | ускорение свободного падения |
| \(\nu\) | \(1.5 \cdot 10^{-5}\ m^2/s\) | кинематическая вязкость воздуха |
| \(deg2rad\) | \(\pi / 180\) | перевод градусов в радианы |
| \(rad2deg\) | \(180 / \pi\) | перевод радиан в градусы |
3. Геометрия крыла и оперения
Нос модели является нулевой точкой конструктора. Координаты деталей \(x_{station}\)
задаются от носа назад, а в физике положение панели относительно ЦТ переводится как:
\[
r_x = x_{CG} - x_{station}
\]
Площадь трапеции
\[
S = b \cdot {c_r + c_t \over 2}
\]
Здесь \(b\) - полный размах, \(c_r\) - корневая хорда, \(c_t\) - концевая хорда.
Удлинение
\[
AR = {b^2 \over S}
\]
Сужение и средняя аэродинамическая хорда
\[
\lambda = {c_t \over c_r}
\]
\[
\bar c = {2 \over 3} c_r {1 + \lambda + \lambda^2 \over 1 + \lambda}
\]
Положение САХ
\[
y_{\bar c} = {b \over 2} {1 + 2\lambda \over 3(1 + \lambda)}
\]
\[
x_{LE,\bar c} = x_{LE} + \tan(\Lambda) y_{\bar c}
\]
\[
x_{AC} = x_{LE,\bar c} + 0.25\bar c
\]
Панелизация
Поверхность делится по полуразмаху на сегменты. Для локальной координаты
\(\eta \in [0,1]\):
\[
c(\eta) = c_r + (c_t - c_r)\eta
\]
\[
\Delta S = \Delta b \cdot c(\eta)
\]
\[
x_{LE,panel} = x - 0.25c_r - \tan(\Lambda){b \over 2}\eta
\]
\[
x_{AC,panel} = x_{LE,panel} + 0.25c(\eta)
\]
4. Центр масс
В расчетном режиме центр масс считается как сумма моментов масс деталей:
\[
x_{CG} = {\sum_i m_i x_i \over \sum_i m_i}
\]
Точки приложения масс:
\[
x_{nose} = noseMassX
\]
\[
x_{fuselage} = 0.48L
\]
\[
x_{wingMass} = x_{LE,\bar c} + 0.45\bar c
\]
\[
x_{tailMass} = x_{LE,\bar c,tail} + 0.45\bar c_{tail}
\]
Если сумма масс деталей больше заданной массы модели, все массы деталей масштабируются:
\[
k = {m_{model} \over \sum_i m_i}
\]
\[
m_i' = k m_i
\]
Если сумма масс деталей меньше массы модели, остаток добавляется как payload в точку носовой массы:
\[
m_{payload} = \max(0, m_{model} - \sum_i m_i)
\]
\[
x_{payload} = noseMassX
\]
В ручном режиме ЦТ задается от передней кромки крыла:
\[
x_{CG} = x_{LE,wing} + x_{CG,fromWingLE}
\]
Положение ЦТ в процентах САХ:
\[
CG_{\%MAC} = 100 {x_{CG} - x_{LE,\bar c} \over \bar c}
\]
5. Статическая устойчивость
Для V-образного хвоста эффективная площадь по тангажу уменьшается через проекцию:
\[
S_{tail,pitch} = S_{tail}\cos^2(\Gamma_V)
\]
Плечо хвоста и хвостовой объем:
\[
l_t = x_{tail} - x_{CG}
\]
\[
V_H = {S_{tail,pitch}\max(l_t, 0) \over S_{wing}\bar c}
\]
Нейтральная точка в текущей упрощенной модели:
\[
x_{NP} = x_{AC,wing} + \bar c \cdot clamp(0.72V_H,\ 0,\ 0.65)
\]
Статический запас:
\[
SM = 100 {x_{NP} - x_{CG} \over \bar c}
\]
6. Воздушная среда
Полная скорость воздуха в точке:
\[
\vec V_{air}(x,y,z,t) = \vec V_{wind} + \vec V_{thermal}(x,y,z) + \vec V_{turb}(x,y,z,t)
\]
\[
\vec V_{wind} = (wind,\ 0,\ 0)
\]
Термики
Для каждого термика с центром \((x_0,z_0)\), радиусом \(R\), высотой \(H\)
и силой \(w_0\):
\[
d = \sqrt{(x - x_0)^2 + (z - z_0)^2}
\]
\[
V_{thermal,y} =
\begin{cases}
w_0 \cdot scale \cdot \left(1 - {d \over R}\right)^2
\left(1 - {y \over H}\right), & d < R,\ 0 \le y < H \\
0, & otherwise
\end{cases}
\]
Турбулентность
\[
h_f = 0.35 + 0.65e^{-\max(y,0)/10}
\]
\[
A = 0.16 \cdot turb \cdot h_f
\]
\[
p_1 = 0.31x + 0.11y + 0.23z + 0.85t
\]
\[
p_2 = 0.17x - 0.19y + 0.29z + 1.20t
\]
\[
p_3 = 0.27x + 0.07y - 0.21z + 0.70t
\]
\[
\vec V_{turb} =
\left(
A\sin p_1 + 0.35A\sin(0.63p_2 + 1.7),\
0.25A\sin(p_2 + 2.1),\
0.60A\sin(p_3 - 0.8)
\right)
\]
Турбулентность задается как скорость воздушной массы, а не как добавочная сила.
Поэтому она не должна работать как скрытая подъемная сила.
7. Кинематика и локальный угол атаки
Начальное состояние
\[
v_x(0) = v_0\cos(\theta_0)
\]
\[
v_y(0) = v_0\sin(\theta_0)
\]
\[
v_z(0) = 0
\]
Скорость точки панели
\[
\vec V_{point} = \vec V_{CG} + \vec \omega_{world} \times \vec r_{world}
\]
\[
\vec V_{rel,world} = \vec V_{point} - \vec V_{air}
\]
\[
\vec V_{rel,body} = q^{-1}\vec V_{rel,world}q
\]
Базис поверхности
Для горизонтальной поверхности до диэдра:
\[
\hat f = (\cos i,\ \sin i,\ 0)
\]
\[
\hat n = (-\sin i,\ \cos i,\ 0)
\]
\[
\hat s = (0,\ 0,\ side)
\]
Для вертикальной поверхности:
\[
\hat f = (\cos i,\ 0,\ \sin i)
\]
\[
\hat n = (-\sin i,\ 0,\ \cos i)
\]
\[
\hat s = (0,\ 1,\ 0)
\]
Угол атаки панели
\[
V_f = \vec V_{rel,body} \cdot \hat f
\]
\[
V_n = \vec V_{rel,body} \cdot \hat n
\]
\[
\alpha = clamp\left(\arctan2(-V_n,\ V_f),\ -60^\circ,\ 60^\circ\right)
\]
8. Аэродинамические коэффициенты
Табличная поляра
Для профиля E387 используется таблица \(C_L\), \(C_D\), \(C_m\).
Между соседними точками применяется линейная интерполяция:
\[
t = {\alpha - \alpha_i \over \alpha_{i+1} - \alpha_i}
\]
\[
C(\alpha) = C_i + (C_{i+1} - C_i)t
\]
За пределами таблицы применяется мягкая экстраполяция:
\[
e = |\alpha_{deg} - \alpha_{edge}|
\]
\[
d = clamp\left(1 - {e \over 45},\ 0.35,\ 1\right)
\]
\[
C_L' = dC_L
\]
\[
C_D' = C_D + 0.0009e^2
\]
\[
C_m' = C_m - 0.0015(\alpha_{deg} - \alpha_{edge})
\]
Если табличный \(C_D\) не содержит индуктивное сопротивление:
\[
C_D = C_{D,profile} + {C_L^2 \over \pi \cdot 0.78 \cdot \max(AR, 0.8)}
\]
Аналитическая модель для простых профилей
\[
\alpha_e = \alpha - \alpha_{0L}
\]
\[
C_{L,linear} = C_{L\alpha}\alpha_e
\]
\[
C_{L,stall} = clamp(C_{L\alpha}\alpha_{stall},\ -C_{L,max},\ C_{L,max})
\]
До срыва:
\[
C_L = clamp(C_{L,linear},\ -C_{L,max},\ C_{L,max})
\]
После срыва:
\[
excess = |\alpha_e| - \alpha_{stall}
\]
\[
d = clamp\left(1 - {excess \over 55^\circ},\ 0.42,\ 1\right)
\]
\[
C_L = sign(\alpha_e)|C_{L,stall}|d
\]
Сопротивление и момент:
\[
C_{D,i} = {C_L^2 \over \pi \cdot 0.75 \cdot \max(AR, 0.8)}
\]
\[
s_d = {\max(0,\ |\alpha_e| - \alpha_{stall}) \over 35^\circ}
\]
\[
C_D = C_{D0} + C_{D,i} + 0.07s_d^2
\]
\[
C_m = C_{m0} - 0.025\alpha_e
\]
9. Силы и моменты
Динамическое давление и число Рейнольдса
\[
q = {1 \over 2}\rho V^2
\]
\[
Re = {Vc \over \nu}
\]
Подъемная сила и сопротивление панели
\[
L = qSC_L
\]
\[
D = qSC_D
\]
Направление подъемной силы берется как нормаль поверхности, очищенная от компоненты вдоль потока:
\[
\hat v = {\vec V_{rel,body} \over |\vec V_{rel,body}|}
\]
\[
\hat l = normalize(\hat n - \hat v(\hat n \cdot \hat v))
\]
Сила в системе тела:
\[
\vec F_{body} = L\hat l - D\hat v
\]
\[
\vec F_{world} = q_{rot}\vec F_{body}q_{rot}^{-1}
\]
Момент панели
\[
\vec M_{Cm,body} = \hat a \cdot qScC_m
\]
\[
\hat a =
\begin{cases}
(0,\ 0,\ 1), & horizontal \\
(0,\ 1,\ 0), & vertical
\end{cases}
\]
\[
\vec \tau_{surface} = \vec r \times \vec F_{world} + \vec M_{Cm,world}
\]
Суммарные сила и момент
\[
\vec F_{total} = (0,\ -mg,\ 0) + \sum_j \vec F_j
\]
\[
\vec \tau_{total} = \sum_j \vec \tau_j
\]
Демпфирование тангажа
\[
M_q = q_{ref}S_{wing}\bar c(-0.08\omega_z)
\]
\[
q_{ref} = {1 \over 2}\rho V_{CG,rel}^2
\]
При отсутствии турбулентности боковые силы дополнительно приглушаются:
\[
k_{lat} =
\begin{cases}
0.45, & turb > 0 \\
0.05, & turb = 0
\end{cases}
\]
\[
F_z \leftarrow k_{lat}F_z,\quad
\tau_x \leftarrow k_{lat}\tau_x,\quad
\tau_y \leftarrow k_{lat}\tau_y
\]
10. Интегратор движения
Шаг симуляции дробится до \(dt \le 0.001s\). Дальше используется явная пошаговая схема.
\[
\vec a = {\vec F_{total} \over m}
\]
\[
\vec v_{t+dt} = \vec v_t + \vec a dt
\]
\[
\vec x_{t+dt} = \vec x_t + \vec v_{t+dt}dt
\]
Тензор инерции в диагональном приближении
\[
h = \max(0.18\bar c,\ 0.015)
\]
\[
I_x = {m(b^2 + h^2) \over 12}
\]
\[
I_y = {m(L^2 + h^2) \over 12}
\]
\[
I_z = {m(L^2 + b^2) \over 12}
\]
Угловое ускорение
\[
\dot\omega_x = 0.22 \cdot clamp\left({\tau_x \over I_x},\ -60,\ 60\right)
\]
\[
\dot\omega_y = 0.22 \cdot clamp\left({\tau_y \over I_y},\ -60,\ 60\right)
\]
\[
\dot\omega_z = clamp\left({\tau_z \over I_z},\ -60,\ 60\right)
\]
Обновление и демпфирование угловых скоростей:
\[
\omega_i \leftarrow clamp(\omega_i + \dot\omega_i dt,\ -8,\ 8)
\]
\[
\omega_x,\omega_y \leftarrow e^{-8dt}\omega_x,\ e^{-8dt}\omega_y
\]
\[
\omega_z \leftarrow e^{-4dt}\omega_z
\]
Кватернион ориентации
\[
\dot q = q \otimes (0,\omega_x,\omega_y,\omega_z)
\]
\[
q_{t+dt} = normalize(q_t + 0.5\dot q dt)
\]
Угол траектории:
\[
\gamma = \arctan2(v_y,\ \sqrt{v_x^2 + v_z^2})
\]
11. Диагностические метрики
12. Оптимизация ЦТ и хвоста
Кнопка оптимизации перебирает ручной \(x_{CG,fromWingLE}\) и угол хвоста.
Целевая функция сейчас - максимальная дальность при фильтре невалидных траекторий.
\[
score =
\begin{cases}
x_{landing}, & valid \\
-\infty, & invalid
\end{cases}
\]
Траектория считается невалидной, если:
\[
\min(v_x) < 0.2
\]
\[
y_{max} > h_0 + 12
\]
\[
pitch_{max} - pitch_{min} > 140^\circ
\]
Поиск идет в два прохода: грубая сетка \(4mm / 0.5^\circ\),
затем точная сетка вокруг лучшего результата \(1mm / 0.1^\circ\).
13. Текущие ограничения
- Нет полноценной нестационарной аэродинамики и динамического срыва.
- Нет гибкости крыла, фюзеляжа и V-оперения.
- Поляра E387 задана для одной области Re, интерполяции по Re пока нет.
- Киль и V-tail сведены к простым панелям без точной интерференции.
- Интегратор простой и малошаговый; это удобно для отладки, но не финальная FDM-схема.
- BOO-пресет пока не валидирован по реальным логам полета.