Разработка математической модели парашютной посадки беспилотного летательного аппарата

1.  Введение

 В данном руководстве представлен обзор поэтапной разработки и верификации модели парашютной посадки беспилотного летательного аппарата (БЛА) самолётной схемы в REPEAT VISION при помощи блока ОДУ, рис. 1:

Рис. 1 Блок ОДУ в REPEAT VISION

 Математическая модель парашютной посадки БЛА основана на работах-монографиях о динамике полёта и управлении БЛА Моисеева В.С. [1],[2] (стр. 176 – 189 в [1]). Верификация модели будет проведена относительно результатов моделирования [1], стр.187.

 Большая часть теоретического материала будет изложена в сокращённом виде, за более подробным описанием математической модели обратитесь к источнику [1],[2].

2.  Описание процесса парашютной посадки БЛА

 Моделирование процесса парашютной посадки БЛА позволяет определить точку его приземления, что сократит время его поиска для повторного применения.

 В типовую моделируемую системы парашютной посадки (СПП) БЛА входят следующие элементы:

• вытяжной парашют (ВП);

• основной парашют (ОП);

• метательный (пиротехнический) заряд (МЗ);

• контейнер для размещения ВП и ОП.

 Процесс парашютной посадки БЛА разбивается на следующие основные этапы:

flowchart TD A(__1__ Выброс ВП с помощью МЗ из контейнера СПП после получения соответствующей команды от САУ, выход в поток и наполнение ВП / __1__ Discharge of EP using PC_2 from the PLS container after receiving the corresponding command from ACS, exit to the flow and filling of EP) --> B("`__2__ Вытягивание купола и строп ОП на всю длину / __2__ Pulling the canopy and MP lines to their full length`") B("`__2__ Вытягивание купола и строп ОП на всю длину / __2__ Pulling the canopy and MP lines to their full length`") --> C(__3__ Наполнение купола ОП до площади ограниченной его рифовкой / <br>__3__ Filling the MP dome to the area limited by its reefing) C(__3__ Наполнение купола ОП до площади ограниченной его рифовкой / <br>__3__ Filling the MP dome to the area limited by its reefing) --> D(__4__ Снижение БЛА на зарифованном куполе ОП / <br>__4__ UAV descent on a reefed MP dome) D(__4__ Снижение БЛА на зарифованном куполе ОП / <br>__4__ UAV descent on a reefed MP dome) --> Е("__5__ Разрифовка и наполнение до полной площади купола ОП / <br>__5__ Reefing and filling to the full area of ​​the MP dome") Е("__5__ Разрифовка и наполнение до полной площади купола ОП / <br>__5__ Reefing and filling to the full area of ​​the MP dome") --> F(__6__ Снижение БЛА на полном куполе ОП, приземление и отцепка ОП от БЛА / <br>__6__ UAV descent with the main canopy fully deployed, landing and detachment of the MP from the UAV)

 Моделирование парашютной посадки БЛА выполняется при следующих допущениях:

• Не учитываются массы заряда и контейнера;

• Применяемые в СПП парашюты являются безынерционными объектами;

• Все действующие в процессе посадки БЛА силы приведены к его центру масс (ЦМ);

• Движение БЛА в процессе посадки рассматривается как движение его ЦМ в вертикальной плоскости;

• Сила лобового сопротивления БЛА нулевая.

3.  Математическая модель процесса парашютной посадки БЛА

 Общие уравнения движения ЦМ БЛА в процессе посадки имеют вид (1) :

(1)\[\begin{split}m\cdot \dot{V} = -G ⋅ sin{\theta} - X - X_{п}\\m\cdot V \cdot\dot{\theta}= -G ⋅ cos{\theta}\\\dot{x}=V\cdot cos{\theta}\\\dot{y}=V\cdot sin{\theta} \end{split}\]

где:

 \(X\) - сила лобового сопротивления БЛА;

 \(X_{п} = \Phi_{вп} \cdot \frac{\rho \cdot V_{1}^2}{2\cdot m}\) - сила лобового сопротивления парашюта;

 \(G\) - сила тяжести БЛА: \(G=m \cdot g\)

 \(V\) - скорость посадки БЛА;

 \(m\) - масса БЛА;

 \(\theta\) - угол посадки БЛА;

 \(x\) - координата БЛА по Х координате;

 \(y\) - координата БЛА по Y координате.

 Так как отсутствуют данные по коэффициенту лобового сопротивления БЛА, \(X\) = 0.

 Начальные условия для этой системы дифференциальных уравнений записываются следующим образом (2) :

(2)\[\begin{split}V(t_{0})=V_{пос}\\\theta(t_{0})=\theta_{пос}\\x(t_{0})=x_{пос}\\y(t_{0})=y_{пос}\end{split}\]

 Достаточно подробный анализ модели приведен в монографии [2].

 При моделировании посадки БЛА рассматриваются следующие основные этапы:

flowchart LR A(__1__ Ввод в действие вытяжного парашюта и вытягивание им основного парашюта / <br>__1__ Deploying the pilot chute and pulling it out of the main parachute) --> B("`__2__ Наполнение основного парашюта / <br>__2__ Inflating the main parachute`") B("`__2__ Наполнение основного парашюта / <br>__2__ Inflating the main parachute`") --> C(__3__ Спуск и приземление БЛА на основном парашюте / <br>__3__ Descent and landing of the UAV using the main parachute)

 Вследствие этого общую систему сведем к трем системам дифференциальных уравнений с взаимосвязанными начальными условиями, каждая из которых описывает соответствующий этап посадки БЛА.

 Модель процесса посадки БЛА на этапе 1 имеет вид (3) :

(3)\[\begin{split}\dot{V_{1}}= -g \cdot sin{\theta_{1}}-\Phi_{вп} \cdot \frac{\rho \cdot V_{1}^2}{2 \cdot m}\\\dot{\theta_{1}}=\frac{-g\cdot cos{\theta_{1}}}{V_{1}}\\\dot{x_{1}}=V_{1}\cdot cos{\theta_{1}}\\\dot{y_{1}}=V_{1}\cdot sin{\theta_{1}}\end{split}\]

где \(\Phi_{вп}\) – эффективное сопротивление ВП.

 Эффективное сопротивление ВП \(\Phi_{вп}\) имеет формулу (4) :

(4)\[\Phi_{вп} =c_{вп}\cdot F_{вп} \]

где:

 \(c_{п}\) – коэффициент сопротивления вытяжного парашюта;

 \(F_{п}\) – площадь купола вытяжного парашюта.

 Начальные условия дифференциальных уравнений первого этапа (5) :

(5)\[\begin{split}V_{1}(t_{0})=V_{пос}\\\theta_{1}(t_{0})=\theta_{пос}\\x_{1}(t_{0})=x_{пос}\\y_{1}(t_{0})=y_{пос}\end{split}\]

 В работах [1], [2] в формуле угловой скорости \(\dot{\theta}\) в правой части уравнения выражение \(G ⋅ cos{\theta}\) без знака минуса, однако при начальном условии \(\theta\) = -180\(^о\) в правой части уравнения будет возникать отрицательная производная, а не положительная, что ведет к убыванию угла, поэтому был добавлен отрицательный множитель.

 Для определения момента времени \(t_{1}\) окончания этапа 1 опишем процесс вытягивания ОП общей массой \(m_{ОП}\) и длиной \(l_{ОП}\) с помощью силы \(X_{ВП}(t)\) дифференциальным уравнением второго порядка (6) :

(6)\[m_{оп}\cdot \ddot{L}=X_{вп}(t) \]

 Начальные условия для этого уравнения имеют вид (7) :

(7)\[\begin{split}L(t_{0})=L_{0}\\\dot{L(t_{0})}=0\end{split}\]

 В уравнении использованы следующие обозначения:

• \(L = L(t)\) – длина вытянутой в момент времени \(t\) части ОП;

• \(m_{ОП}= m_{ОП}(L)\) – масса вытянутой в этот момент времени части ОП длиной \(L\).

 Массу ОП, вытянутого на длину \(L\), можно приближенно вычислить по формуле (8) :

(8)\[m_{ОП}(L)=\frac{\pi \cdot \mu \cdot d^2 \cdot L}{4}\]

где:

\(\mu\) – средняя плотность материала ОП;

\(d\) – средний диаметр ОП.

 Используя выражение (9) :

(9)\[X_{вп}(t)=\frac{\Phi_{вп} \cdot \rho \cdot V_{1}^2}{2} \]

и формулу массы ОП ДУ второго порядка (5) можно представить виде системы двух уравнений первого порядка (10) :

(10)\[\begin{split}\dot{L}=l\\\dot{l}=\frac{2 \cdot \Phi_{вп} \cdot \rho \cdot V_{1}^2}{\pi \cdot \mu \cdot d^2 \cdot L}\end{split}\]

с начальными условиями (11) :

(11)\[\begin{split}L(t_{0})=L_{0}\\l(t_{0})=0\end{split}\]

 Тогда момент времени \(t_{1}\) будет определяться из условия (12) :

(12)\[L(t_{1}) = L_{оп}\]

т.е. когда вытяжной парашют вытягивается на длину основного парашюта, раскрывая его.

 Этап 2 посадки БЛА начинается при раскрытии основного парашюта и описывается следующей системой дифференциальных уравнений (13) :

(13)\[\begin{split}\dot{V_{2}}= -g \cdot sin{\theta_{2}}-\Phi_{оп}(t,V) \cdot \frac{\rho \cdot V_{2}^2}{2 \cdot m}\\\dot{\theta_{2}}=\frac{-g\cdot cos{\theta_{2}}}{V_{2}}\\\dot{x_{2}}=V_{2} \cdot \cos{\theta_{2}}\\\dot{y_{2}}=V_{2}\cdot sin{\theta_{2}}\end{split}\]

с начальными условиями (14) :

(14)\[\begin{split}V_{2}(t_{1})=V_{1}(t_{1})\\\theta_{2}(t_{1})=\theta_{1}(t_{1})\\x_{2}(t_{1})=x_{1}(t_{1})\\y_{2}(t_{1})=y_{1}(t_{1})\end{split}\]

\(\Phi_{оп}(t,V)\) вычисляется по формуле (15) :

(15)\[\Phi_{оп}(t,V)= \Phi_{max} \cdot (\frac{t\cdot V}{n\cdot D_{оп}})^k\]

где:

 \(\Phi_{max} = c_{оп} \cdot F_{оп}\) - максимальное значение эффективной площади основного парашюта;

 \(D_{оп}\) - номинальный диаметр парашюта;

 \(k, n\) - эмпирические коэффициенты.

 Этап 3 предлагается описывать моделью вида (16) :

(16)\[\begin{split}\dot{V_{3}}= -g \cdot sin{\theta_{3}}-\Phi_{max} \cdot \frac{\rho \cdot V_{3}^2}{2 \cdot m}\\\dot{\theta_{3}}=\frac{-g\cdot cos{\theta_{3}}}{V_{3}}\\\dot{x_{3}}=V_{3}\cdot cos{\theta_{3}}\\\dot{y_{3}}=V_{3}\cdot sin{\theta_{3}}\end{split}\]

с начальными условиями (17) :

(17)\[\begin{split}V_{3}(t_{2})=V_{2}(t_{2})\\\theta_{3}(t_{2})=\theta_{2}(t_{2})\\x_{3}(t_{2})=x_{2}(t_{2})\\y_{3}(t_{2})=y_{2}(t_{2})\end{split}\]

 Начинается третий этап при \(\Phi_{оп}(t,V)=\Phi_{max}\), т.е. когда множитель \(u=\frac{t \cdot V}{n \cdot D_{оп}}=1\).

4.  Исходные данные

 Пусть БЛА, взлетевший с МПУ после выполнения полетного задания, осуществляет парашютную посадку в стартовой СК, представленной на рис. 2:

Рис. 2 Схема парашютной посадки после полетного задания

 Считается, что посадка БЛА происходит при нормальных атмосферных условиях (\(р_{0}\) = 760 мм.рт.ст., \(\tau_{0}\) = +15\(^о\)C) и отсутствии ветра.

 Исходные данные:

• Беспилотный летательный аппарат:

Наименование характеристики	Значение
Посадочная масса m	300 [кг]

• Вытяжной парашют (ВП):

Наименование характеристики	Значение
Площадь купола F_{вп}	1.02 [м]
Коэффициент сопротивления c_{вп}	0.86

• Основной парашют (ОП):

Наименование характеристики	Значение
Номинальный диаметр D_{оп}	14 [м]
Площадь купола F_{оп}	312.4 [м]
Коэффициент сопротивления c_{оп}	0.81
Эмпирический коэффициент k	1.6
Эмпирический коэффициент n	2

• Укладка ОП:

Наименование характеристики	Значение
Длина вытянутого ОП L_{оп}	31 [м]
Средний диаметр укладки d_{оп}	0.25 [м]
Удельный вес укладки из парашютного шелка μ	1.36

• Остальные параметры посадки:

Наименование характеристики	Значение
Плотность воздуха ρ	1.225 кг/м\(^3\)
Ускорение свободного падения g	9.81 м/с\(^2\)

 Начальные условия ОДУ (1 этап) (18) :

(18)\[\begin{split}V_{1}(0)=50 \frac{м}{с}\\θ_1 (0)=-180^о\\x_{1}(0)=1000 \hspace{1mm} м\\y_{1}(0)=800\\L(0) = 0.2\hspace{1mm}м\\l = 0\end{split}\]

5.  Разработка модели процесса посадки БЛА в REPEAT VISION

 В данном разделе представлен обзор поэтапной разработки рассматриваемой математической модели при помощи блока ОДУ.

Подсказка

Также готовая модель представлена в демо проекте «[DEMO] Парашютная посадка беспилотного летательного аппарата» в разделе Демо проекты.

 Блок ОДУ [3] предназначен для численного решения систем обыкновенных дифференциальных уравнений (ОДУ), выраженных в форме Коши (относительно старшей производной): \(y'' = f(t,y,y')\). При помощи этого блока можно описать комплексную динамическую систему.

 Весь процесс настройки блока ОДУ состоит из 5 ключевых шагов:

flowchart LR A(__1__ Очистка шаблона и подготовка рабочего пространства в редакторе кода) --> B("`__2__ Добавление входных/выходных портов и параметров`") B("`__2__ Добавление входных/выходных портов и параметров`") --> C(__3__ Инициализация констант и начальных условий ОДУ) C(__3__ Инициализация констант и начальных условий ОДУ) --> D(__4__ Написание системы ОДУ) D(__4__ Написание системы ОДУ) --> Е("__5__ Настройка и модификация проекта")

 Поэтапно опишем процесс работы с блоком:

1. При создании и нажатии на блок ОДУ появляется редактор скрипта по написанию системы ОДУ, рис. 3:

Рис. 3 Код шаблонного блока ОДУ

 В шаблонном блоке ОДУ представлена модель отскока мяча от двух стенок, для написания нашей модели очистим полностью от кода разделы Инициализация, Основной блок, Дополнительные блоки, рис. 4:

Рис. 4 Редактор кода блока ОДУ после очистки

далее удалим переменные, находящиеся в правой панели (необходимо навести на параметр и нажать значок урны), рис. 5:

Рис. 5 Удаление параметров шаблона блока ОДУ

удалим в блоке ОДУ в разделе (1) существующие события, раскрыв вкладку существующего события и нажав на кнопку Удалить (2), рис. 6:

Рис. 6 Удаление событий шаблона блока ОДУ

 Блок ОДУ теперь готов к описанию математической модели динамической системы парашютной посадки БЛА.

2. Перед написанием системы ОДУ необходимо в блоке добавить:

• параметры входных и выходных портов;

• параметры, участвующие в процессе расчёта ОДУ;

• параметры интегрирования.

 Параметры входного порта:

1. timer – для снятия значения времени моделирования, участвует в формуле \(\Phi_{оп}(t,V)\).

 Параметры выходных портов (см. рис. 7 ниже):

1. velocity_sensor – для снятия значения параметра скорости посадки БЛА;

2. theta_sensor – для снятия значения параметра угла БЛА, [рад];

3. x_sensor – для снятия значения координаты x, [м];

4. y_sensor – для снятия значения координаты y, [м];

5. L_OP_sensor – для снятия значения длины вытянутой части ОП, [м].

Рис. 7 Требуемые параметры портов блока ОДУ

 Параметры (см. рис. 8 ниже):

 m (double) – масса БЛА;

 c_VP (double) – коэффициент сопротивления вытяжного парашюта;

 f_VP (double) – площадь купола вытяжного парашюта;

 ro (double) – плотность воздуха;

 mu (double) – удельный вес укладки из парашютного шелка;

 k (double) – эмпирический коэффициент при вычислении \(\Phi_{оп}(t,V)\);

 n (double) – эмпирический коэффициент при вычислении \(\Phi_{оп}(t,V)\);

 g (double) – ускорение свободного падения;

 phi_VP (double) – эффективное сопротивление ВП;

 phi_OP (double) – эффективное сопротивление ОП;

 u (double) – множитель при расчете \(\Phi_{оп}\);

 gndHit (int) – переменная-сигнал при приземлении БЛА;

 c_OP (double) – коэффициент сопротивления основного парашюта;

 f_OP (double) – площадь купола основного парашюта;

 D_OP (double) – номинальный диаметр ОП;

 d_OP (double) – средний диаметр укладки ОП.

Рис. 8 Требуемые параметры блока ОДУ

 Параметры интегрирования:

 velocity – основной параметр для интегрирования производной ускорения посадки БЛА \(\dot{V}\), необходимый для вычисления скорости \(V\);

 theta – параметр для интегрирования производной угла \(\dot{\theta}\);

 x – параметр для интегрирования производной координаты \(\dot{x}\);

 y – параметр для интегрирования производной координаты \(\dot{y}\);

 l_OP – параметр для интегрирования производной ускорения раскрытия основного парашюта;

 L_OP – параметр для интегрирования производной скорости раскрытия основного парашюта для определения ее длины \(L_{оп}\).

 Интегрируемые параметры - состояние разработанной системы. Блок ОДУ будет интегрировать их производные, чтобы получить численное решение системы ОДУ, рис. 9:

Рис. 9 Параметры интегрирования блока ОДУ

3. После добавления параметров, в редакторе кода переходим в блок Инициализация, где созданным постоянным параметрам присваиваем значения и формируем начальные условия системы ОДУ, т.е. задаем начальные значения интегрируемых параметров.

m = 300
c_VP = 0.86
f_VP = 1.02
phi_VP = c_VP*f_VP
d_OP = 0.25
D_OP = 14
f_OP = 312.4
c_OP = 0.81
mu = 1.36
k = 1.6
n = 2
g = 9.81
ro = 1.225
gndHit = 0
velocity = 50
theta = -Math.PI
x = 1000
y = 800
L_OP = 0.2
l_OP = 0

• Строки 1-14. Присваиваются значения параметрам блока ОДУ.

• Строки 15-20. Инициализация интегрируемых параметров: начальная скорость посадки БЛА (15), начальный угол при посадке БЛА (16), начальная координата по X (17), начальная координата по Y (18), начальная длина ОП (19)

4. Переходим к основной части разработки математической модели посадки БЛА: необходимо описать систему ОДУ, реализуя функции \(f(t,y)\), в основном блоке. Начнём с первого этапа посадки БЛА, когда дифференциальные уравнения описывают процесс работы вытяжного парашюта:

velocity' = -g*Math.sin(theta) -phi_VP*velocity*velocity*ro/2/m;
theta' = -g*Math.cos(theta)/velocity;
x'=velocity*Math.cos(theta);
y'=velocity*Math.sin(theta);
l_OP'=2*phi_VP*velocity*velocity*ro/(Math.PI*mu*d_OP*d_OP*L_OP);
L_OP'= l_OP;

 Второй этап возникает при полном раскрытии основного парашюта (\(L_{оп} \geq31\)), когда силу сопротивления падению создает ОП:

velocity' = -g*Math.sin(theta) - phi_OP*velocity*velocity*ro/2/m ;
theta' = -g*Math.cos(theta)/velocity;
x'=velocity*Math.cos(theta);
y'=velocity*Math.sin(theta);
l_OP'= 0;
L_OP'= 0;

 Блок ОДУ поддерживает работу условных операторов if-else, воспользуемся данной функцией для написания динамической системы. Дифференциальные уравнения угла и координат не изменяются, их можно вынести в общую область кода. Код в основном блоке будет выглядеть как:

if(L_OP_sensor<31){
velocity' = -g*Math.sin(theta) -phi_VP*velocity*velocity*ro/2/m ;
l_OP'=2*phi_VP*velocity*velocity*ro/(Math.PI*mu*d_OP*d_OP*L_OP);
L_OP'= l_OP;
} else {
velocity' = -g*Math.sin(theta) - phi_OP*velocity*velocity*ro/2/m ;
l_OP'= 0;
L_OP'= 0;
}
theta' = -g*Math.cos(theta)/velocity;
x'=velocity*Math.cos(theta);
y'=velocity*Math.sin(theta);

• Строки 1-5. Блок кода, с уравнениями первого этапа посадки, когда работает вытяжной парашют. Данный блок выполняется, пока не раскроется основной парашют на всю длину (L_OP_sensor<31);

• Строки 5-9. Блок кода, с уравнениями второго этапа посадки при раскрытии основного парашюта. Расчет модели в этом блоке происходит, когда основной парашют вытянут на всю длину (L_{оп} = 31);

• Строка 10. Уравнение скорости изменения угла посадки;

• Строка 11. Уравнение скорости изменения координаты X;

• Строка 12. Уравнение скорости изменения координаты Y.

 Однако код расчета \(\Phi_{оп}\) еще не был продемонстрирован, в том числе и передача значений проинтегрированных параметров в выходные порты: всё это происходит во вкладке Дополнительно в подразделе После интегрирования. Здесь также происходит расчет 3 этапа, когда эффективная площадь сопротивления основного парашюта \(\Phi_{оп}\) достигает максимального значения \(\Phi_{макс.}= c_{оп} \cdot F_{оп}\) при \(u=1\).

velocity_sensor = velocity;
theta_sensor = theta;
x_sensor = x;
y_sensor = y;
L_OP_sensor = L_OP;
if(L_OP>31){
L_OP = 31;
L_OP_sensor = L_OP;
}
u = (timer*velocity)/(n*D_OP);
if(u>1){ 
u = 1;
}
phi_OP = (c_OP*f_OP) * Math.pow(u, k);

• Строки 1-5. Вывод значений интегрируемых параметров в выходные порты;

• Строки 6-9. Блок кода для ограничения длины ОП;

• Строка 10. Расчет вспомогательной функции \(u(t,V)\) для расчета силы сопротивления основного парашюта;

• Строки 11-13. Код блока для ограничения \(u(t,V)\);

• Строка 14. Расчет силы сопротивления основного парашюта.

 Параметры интегрирования можно оставить шаблонные, т.е.:

• Максимальное кол-во обращений к функции – 10000;

• Тип интегратора – DormandPrince853;

• Минимальный шаг интегратора – 1e-34;

• Максимальный шаг интегратора – 1;

• Абсолютная точность – 1e-6;

• Относительная точность – 1e-6.

 Код в редакторе блока ОДУ:

Инициализация

m=300;
c_VP = 0.86;
f_VP = 1.02;
phi_VP =c_VP*f_VP;
d_OP = 0.25;
D_OP = 14;
f_OP = 312.4;
c_OP = 0.81;
mu = 1.36;
k = 1.6;
n = 2;
g = 9.81;
ro = 1.225;
gndHit = 0;
L_OP = 0.2;
l_OP = 0; 
velocity = 50;
theta = -Math.PI;
x = 1000;
y = 800;

Основной блок

if(L_OP_sensor<31){
velocity' = -g*Math.sin(theta) -phi_VP*velocity*velocity*ro/2/m ;
l_OP'=2*phi_VP*velocity*velocity*ro/(Math.PI*mu*d_OP*d_OP*L_OP);
L_OP'= l_OP;
} else {
velocity' = -g*Math.sin(theta) - phi_OP*velocity*velocity*ro/2/m ;
l_OP'= 0;
L_OP'= 0;
}
theta' = -g*Math.cos(theta)/velocity;
x'=velocity*Math.cos(theta);
y'=velocity*Math.sin(theta);

Дополнительные блоки – Перед интегрированием

- (Пустое поле)

Дополнительные блоки – После интегрирования

velocity_sensor = velocity;
theta_sensor = theta;
x_sensor = x;
y_sensor = y;
L_OP_sensor = L_OP;
if(L_OP_sensor>31){
L_OP = 31;
L_OP_sensor = L_OP;
}
u = (timer*velocity)/(n*D_OP);
if(u>1) { 
u = 1; 
}
phi_OP = (c_OP*f_OP)*Math.pow(u, k);

5. Осталось модифицировать сам проект, добавив блоки:

• Время моделирования (в разделе Симуляция), который подключается ко входному порту блока ОДУ и передает значения времени расчета модели;

• Единицы измерения угла (в библиотеке Автоматика, подраздел Функции) для перевода значения угла посадки БЛА из радиан в градусы, рис. 10:

Рис. 10 Проект с математической моделью посадки БЛА

 Время моделирования проекта выберем равным 200 [с], шаг интегрирования и дискретизации 1 [мс], рис. 11:

Рис. 11 Параметры расчета проекта

 Модель готова к расчету. После получения результатов моделирования можно приступить к верификации расчетов.

6.  Результаты моделирования, верификация расчетов

Таблица 1 – Сравнительный анализ параметров при спуске БЛА

Этап 1

t, с	Скорость V, м/с Эталон	Скорость V, м/с REPEAT VISION	Скорость V, м/с Абс. погр.	Угол \(\theta\), град Эталон	Угол \(\theta\), град REPEAT VISION	Угол \(\theta\), град Абс. погр.	Координата y, м Эталон	Координата y, м REPEAT VISION	Координата y, м Абс. погр.
0.005	49.976	49.978	0.002	-179.94	-179.944	0.00379	800	800	0
0.011	49.953	49.951	0.002	-179.881	-179.876	0.00472	799.999	799.999	0
0.016	49.929	49.929	0.000	-179.821	-179.820	0.001	799.999	799.999	0
0.021	49.906	49.907	0.001	-179.761	-179.764	0.0027	799.998	799.998	0
0.027	49.882	49.880	0.002	-179.702	-179.696	0.00588	799.997	799.996	0.001
0.064	49.72	49.719	0.001	-179.283	-179.279	0.00448	799.98	799.98	0
0.069	49.698	49.698	0.000	-179.223	-179.222	0.001	799.977	799.977	0
0.074	49.675	49.676	0.001	-179.163	-179.165	0.00244	799.973	799.973	0
0.08	49.652	49.651	0.002	-179.103	-179.098	0.00547	799.969	799.969	0
0.496	19.691	19.693	0.002	-171.352	-171.350	0.00187	799.131	799.13	0.001

Этап 2

t, с	Скорость V, м/с Эталон	Скорость V, м/с REPEAT VISION	Скорость V, м/с Абс. погр.	Угол \(\theta\), град Эталон	Угол \(\theta\), град REPEAT VISION	Угол \(\theta\), град Абс. погр.	Координата y, м Эталон	Координата y, м REPEAT VISION	Координата y, м Абс. погр.
1.329	8.904	8.904	0.000	-136.43	-136.404	0.02587	795.256	795.252	0.0041
2.161	7.25	7.252	0.002	-107.391	-107.390	0.00086	789.634	789.632	0.0017
2.994	6.276	6.276	0.000	-95.2	-95.196	0.00405	784.127	784.123	0.00414
3.826	5.59	5.591	0.001	-91.306	-91.306	0.000275	779.212	779.210	0.0019
4.659	5.097	5.098	0.001	-90.282	-90.282	0.00028	774.773	774.769	0.00426
5.492	4.724	4.724	0.000	-90.053	-90.053	0.00018	770.691	770.684	0.0068
6.324	4.427	4.428	0.001	-90.009	-90.009	9E-05	766.885	766.881	0.0038

Этап 3

t, с	Скорость V, м/с Эталон	Скорость V, м/с REPEAT VISION	Скорость V, м/с Абс. погр.	Угол \(\theta\), град Эталон	Угол \(\theta\), град REPEAT VISION	Угол \(\theta\), град Абс. погр.	Координата y, м Эталон	Координата y, м REPEAT VISION	Координата y, м Абс. погр.
18.059	4.358	4.358	0.000	-90	-90	0	715.747	715.730	0.0174
41.528	4.357	4.358	0.001	-90	-90	0	613.497	613.462	0.03515
64.997	4.357	4.358	0.001	-90	-90	0	511.248	511.194	0.05388
88.467	4.357	4.358	0.001	-90	-90	0	408.998	408.922	0.07597
111.936	4.357	4.358	0.001	-90	-90	0	306.749	306.654	0.0947
135.405	4.357	4.358	0.001	-90	-90	0	204.499	204.387	0.11245
158.874	4.357	4.358	0.001	-90	-90	0	102.25	102.119	0.13119
170.609	4.357	4.358	0.001	-90	-90	0	51.125	50.983	0.142235

Таблица 2 – Сравнительный анализ конечных параметров при приземлении БЛА

Параметр	Эталонное решение	REPEAT VISION	Относит. погрешность, %
Время приземления \(t_{пр}\),с	182.34	182.309	0.017
Скорость приземления \(V_{пр}\), м/с	4.357	4.358	0.02
Угол \(\theta_{пр}\), \(^о\)	-90	-90	0
Координата \(x_{пр}\), м	968.5	968.48	0.002

 Результаты моделирования в REPEAT VISION показали высокую точность расчета математической модели посадки БЛА при помощи блока ОДУ, максимальная абсолютная погрешность составила 0.14 [м] (0.28 % в отн. погр.) при расчете высоты БЛА Y.

7.  Модификация модели посадки – столкновение БЛА с землей

 Для демонстрации способности блока ОДУ решать жесткие задачи и добавления визуальной реалистичности, модель может быть дополнена четвертым этапом - упругим столкновением БЛА с земной поверхностью. Данная модификация реализует классическую задачу отскока мяча с энергетическими потерями.

 Столкновение моделируется как пружинно-демпферная система:

• Жесткость пружины (stiffness = 1e4 [Н/м]) - определяет упругость отскока;

• Демпфирование (damping = 50 [Н·с/м]) - обеспечивает постепенную потерю энергии.

 Уравнение движения при контакте с землей (19) :

(19)\[\dot{V} = g + k\cdot y - c \cdot V\]

где:

 \(k\) - коэффициент жесткости;

 \(c\) - коэффициент демпфирования;

 \(y\) - глубина проникновения в землю (координата y).

Важно

Так как момент столкновения (координата \(y<0\)) – это резкое, прерывистое событие, продолжение расчета с прежними производными может привести к нестабильному расчету решателя ОДУ, поэтому необходимо добавить событие (см. рис. 12 ниже), которое при изменении знака координаты \(y\) производило сброс производных. Для этого во вкладке События добавим простое событие (1), раскроем созданную вкладку события и напишем условие с интегрируемым параметром y; (обязательно с ;!) (2) и в выпадающем списке 3 выберем действие Сбросить производное. Написав условие y;, мы даем инструкцию решателю ОДУ, что необходимо произвести сброс значения производных при смене знака высоты \(y\).

Рис. 12 Добавление нового простого события

 Помимо нового события, был модифицирован код для обработки столкновения, изменения коснулись разделов Основной блок и После интегрирования в Дополнительных блоках:

Основной блок

if(L_OP_sensor<31) {
velocity' = -g*Math.sin(theta) -phi_VP*velocity*velocity*ro/2/m ;
l_OP'=2*phi_VP*velocity*velocity*ro/(Math.PI*mu*d_OP*d_OP*L_OP);
L_OP'= l_OP;
theta' = -g*Math.cos(theta)/velocity;
} else {
if(y>0) {
velocity' = -g*Math.sin(theta) - phi_OP*velocity*velocity*ro/2/m ;
theta' = -g*Math.cos(theta)/velocity;
if(gndHit==1){
velocity' = -g*Math.sin(theta);
}
} else {
double stiffness =1e4; 
double damping = 50;
velocity' = g  + stiffness*y - damping*velocity;
theta' = 0;
} 
l_OP'= 0;
L_OP'= 0;
}
x'=velocity*Math.cos(theta);
y'=velocity*Math.sin(theta);

• Строки 1-6. Блок кода, с уравнениями первого этапа посадки, когда работает вытяжной парашют. Данный блок выполняется, пока не раскроется основной парашют на всю длину (L_OP_sensor<31);

• Строки 7-13. Блок кода, с уравнениями второго этапа посадки при раскрытии основного парашюта. Расчет модели в этом блоке происходит, когда основной парашют вытянут на всю длину (L_{оп} = 31). В подблоке 10-13 написано уравнение скорости после столкновения с землей, когда сила сопротивления основного парашюта становится нулевой;

• Строки 14-16. Блок кода с уравнениями четвертого этапа, расчет модели при столкновении с землей.

Дополнительные блоки – После интегрирования

velocity_sensor = velocity;
theta_sensor = theta;
x_sensor = x;
y_sensor = y;
L_OP_sensor = L_OP;
if(L_OP_sensor>31) {
L_OP = 31;
L_OP_sensor = L_OP;
}
u = (timer*velocity)/(n*D_OP);
if(u>1) { 
u = 1;
}
phi_OP = (c_OP*f_OP)*Math.pow(u, k);
if(y_sensor<0){
gndHit =1;
} 

• Строки 15-17. В подразделе После интегрирования добавляются строки 15-17, в которых переменная gndHit сигнализирует о столкновении с землей и становится равной 1.

8.  Результаты моделирования модифицированной модели посадки БЛА

 Графики высоты и скорости БЛА в до и после столкновения с землей имеют вид, рис. 13:

Рис. 13 Графики высоты и скорости посадки БЛА

 Результаты расчета модифицированной модели демонстрируют, что блок ОДУ может работать не только непрерывными процессами, но и с жёсткими системами, где параметры ДУ меняются скачкообразно.

9.  Выводы

 В рамках данного руководства была успешно разработана, реализована и верифицирована математическая модель парашютной посадки беспилотного летательного аппарата самолётной схемы на платформе REPEAT VISION с использованием блока ОДУ.

1. Достигнутые результаты и верификация модели

Проведенное исследование демонстрирует высокую эффективность и точность применения блока ОДУ для решения сложных динамических задач. Результаты расчёта, полученные на платформе REPEAT VISION, с высокой точностью соответствуют данным, приведённым в монографиях (В.С. Моисеев [1] , [2]), служат убедительным доказательством корректности модели.

2. Функциональность и надежность блока ОДУ

В ходе работы была подтверждена способность блока ОДУ платформы REPEAT VISION устойчиво решать системы ОДУ с прерывистой динамикой.

3. Практическая значимость и область применения

Разработанная модель представляет собой готовый инструмент для проектных расчётов, позволяющий с высокой достоверностью определять траекторию и точку приземления БЛА.

10.  Список литературы

---

[1] (1,2,3,4,5,6)

Моисеев В.С. Динамика полета и управление беспилотными летательными аппаратами: монография. – Казань: Редакционно-издательский центр «Школа», 2017. 416 с.

[2] (1,2,3,4,5)

Моисеев В.С. Прикладная теория управления беспилотными летательными аппаратами: монография. Казань: Республиканский центр мониторинга качества образования (РЦ МКО), 2013. 768 с.

[3]

Справка платформы REPEAT VISION: блок ОДУ.