Динамика вертикального полёта летательного аппарата легче воздуха

в 17:21, , рубрики: python, Алгоритмы, вытаскивающая сила, летательный аппарат, математика, разработка под windows, скорость подъёма, уравнение движения

Введение

Определение скорости подъёма и спуска летательных аппаратов легче воздуха (ЛАЛВ) до настоящего времени является практически важной задачей, возникающей при проектировании таких аппаратов.

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

Последнее утверждение требует определённого пояснения, поскольку искушённый читатель хорошо помнит школьный курс физики, в котором решались задачи на высоту подъёма и другие параметры воздушных шаров, заполненных газами легче воздуха или самим подогреваемым во время полёта воздухом.

Все указанные задачи были основаны на равенстве двух сил: силы веса и выталкивающей силы. Газы считались идеальными и их параметры вычислялись по закону Менделеева Клапейрона. Однако, даже простой учёт третьей силы сопротивления воздуха уже приводит к системе дифференциальных уравнений, которая аналитически не решается. Необходимо так же учитывать изменение плотности атмосферного воздуха с высотой подъёма и температурой.

Кроме этого, если нужно рассмотреть не только подъём, но и зависание шара и его спуск на землю, то совсем уж не детская задача получается. Надеюсь, что рассмотрение решения подобной задачи средствами Python не только будет способствовать расширению знаний по физике, но и популяризации самого языка программирования Python. Что я и пытаюсь делать в своих публикациях на этом ресурсе.

Математическая модель полёта ЛАЛВ с оболочкой в форме шара, объём которого не изменяется с изменением высоты

Ограничимся рассмотрением движения его центра масс под действием следующих сил: силы тяжести (G), архимедовой силы (Fa) и силой аэродинамического сопротивления (Fc). Запишем соотношения для определения сил через параметры движения и воздушной среды[3]:

Динамика вертикального полёта летательного аппарата легче воздуха - 1
Динамика вертикального полёта летательного аппарата легче воздуха - 2
Динамика вертикального полёта летательного аппарата легче воздуха - 3

В приведенных формулах приняты обозначения: h – высота подъема шара, dh/dt – вертикальная скорость, m – масса, g – ускорение свободного падения, W – объем шара, c – коэффициент лобового сопротивления, S – характерная площадь сопротивления (площадь Миделя).
Зависимость плотности воздуха от высоты будем полагать экспоненциальной:

Динамика вертикального полёта летательного аппарата легче воздуха - 4

где $ρ_0$ – плотность воздуха на нулевой высоте, b–коэффициент. Сила тяжести направлена вниз, архимедова сила – вверх, а сила аэродинамического сопротивления всегда направлена «против движения», поэтому корректный учет этой силы в уравнениях движения требует введения множителя $-sing(dh/dt)$.

Однако, для наших целей этот факт не имеет принципиального значения, и мы ограничимся рассмотрением только этапа подъема шара, когда сила аэродинамического сопротивления направлена вниз и, следовательно, будет учтена в уравнениях движения со знаком минус. Теперь уравнение движения может быть записано в виде:

Динамика вертикального полёта летательного аппарата легче воздуха - 7, (1)

Дополнительно предположим, что воздушный шар представляет собой однородное тело радиуса R с плотностью $ρ_b$. Тогда величина площади, определяющая его аэродинамическое сопротивление, определится как Динамика вертикального полёта летательного аппарата легче воздуха - 9, объем как Динамика вертикального полёта летательного аппарата легче воздуха - 10, а масса, соответственно, как Динамика вертикального полёта летательного аппарата легче воздуха - 11.

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

Динамика вертикального полёта летательного аппарата легче воздуха - 12, (2)

Введём обозначения:

Динамика вертикального полёта летательного аппарата легче воздуха - 13;
Динамика вертикального полёта летательного аппарата легче воздуха - 14;

и перепишем (2) в виде следующей системы линейных уравнений:

Динамика вертикального полёта летательного аппарата легче воздуха - 15, (3)

Влияние на скорость и высоту подъёма ЛАЛВ температуры атмосферного воздуха

Для этого сначала решим систему (3) с использованием следующего соотношения для зависимости плотности атмосферного воздуха от высоты без учёта температуры:

Динамика вертикального полёта летательного аппарата легче воздуха - 16

Повторим решение системы (3), но уже с использованием соотношения для зависимости плотности воздуха от высоты и температуры:

Динамика вертикального полёта летательного аппарата легче воздуха - 17

где: b=0.000125 — константа, связанная с плотностью воздуха в 1/м.;
a=0.0065 — константа, связанная с температурой воздуха в K/м.
$T_0=300 K$ – температура на уровне моря.

Листинг программы

 # -*- coding: utf8 -*-
from numpy import*
from scipy.integrate import odeint
import matplotlib.pyplot as plt
g=9.81# ускорение свободного падения на земле в м/с2.
rv=1.29# плотность атмосферного воздуха в кг/м3.
rg=0.17# плотность гелия в кг/м3.
R=8# радиус оболочки ЛАЛВ в м.
b=0.000125# константа, связанная с плотностью воздуха в 1/м
a=6.5*10**-3# константа, связанная с температурой воздуха в К/м
c=0.4#коэффициент лобового сопротивления
mo=240#масса в кг
V=(4/3)*pi*R**3
rs=rg+mo/V# суммарная плотность материала ЛАЛВ, массы гелия, и нагрузки
p1=rv/rs# введенный параметр
p2=3*c/(8*R)# введенный параметр
T0=300
def fun(y, t):
         y1, y2= y    
         return [y2,-g+g*p1*exp(-b*y1*T0/(T0-a*y1))-p1*p2*exp(-b*y1*T0/(T0-a*y1))*y2**2]
t =arange(0,1100,0.01)
y0 = [0.0,0.0]
[y1,y2]=odeint(fun, y0,t, full_output=False).T
plt.title("Характеристики  подъёма ЛАЛВ  n Объём: %s м3. Масса : %s кг. n Подъёмная сила: %s kН. "%(round(V,0),mo,round(0.001*g*rv*V,0)))
plt.plot(t/60,y1,label='Максимальная высота подъёма: %s км. n Максимальная скорость: % s м/с .n С учётом температуры воздуха'%(round(max(y1)/1000,2),round(max(y2),2)))
def fun(y, t):
         y1, y2= y
         return [y2,-g+g*p1*exp(-b*y1)-p1*p2*exp(-b*y1)*y2**2]
[y1,y2]=odeint(fun, y0,t, full_output=False).T
plt.plot(t/60,y1,label='Максимальная высота подъёма: %s км. n Максимальная скорость: % s м/с n Без учёта температуры воздуха'%(round(max(y1)/1000,2),round(max(y2),2)))
plt.ylabel('Высота в м')
plt.xlabel(' Время в минутах')
plt.legend(loc='best')
plt.grid(True)
plt.show()

Получим:

Динамика вертикального полёта летательного аппарата легче воздуха - 19

Расчётное значение высоты подъёма ЛАЛВ с учётом температуры меньше, чем без учёта. Скорость подъёма аппарата при этом остаётся неизменной.

Определение характеристик всех фаз полёта ЛАЛВ от старта до приземления

Для построения программы полёта ЛАЛВ рассмотрим условия для следующих периодов времени:

ПодъёмДинамика вертикального полёта летательного аппарата легче воздуха - 20;
Зависание Динамика вертикального полёта летательного аппарата легче воздуха - 21;
ПриземлениеДинамика вертикального полёта летательного аппарата легче воздуха - 22.

Листинг программы

# -*- coding: utf8 -*-
from numpy import*
from scipy.integrate import odeint
import matplotlib.pyplot as plt
g=9.81# ускорение свободного падения на земле в м/с2.
rv=1.29# плотность атмосферного воздуха в кг/м3.
rg=0.17# плотность гелия в кг/м3.
R=8# радиус оболочки стратостата в м.
b=0.000125# константа, связанная с плотностью воздуха в 1/м
a=6.5*10**-3# константа, связанная с температурой воздуха в К/м
c=0.4# коэффициент лобового сопротивления
mo=240# масса в кг
V=(4/3)*pi*R**3
p2=3*c/(8*R)# введенный параметр
T0=300# температура на уровне моря
tz=4000# время зависания в секундах
rgu=1.2# плотность образовавшейся газовой смеси после  стравливания гелия в кг/м3 
tz=4000# время зависания
def fun(y, t):
         y1,y2= y
         if y2<=0:
                  if t<tz:
                            return [y2,-g+g*(rv/(rg+mo/V))*exp(-b*y1*T0/(T0-a*y1))+(rv/(rg+mo/V))*p2*exp(-b*y1*T0/(T0-a*y1))*y2**2]
                  elif t>=tz:
                           return [y2,-g+g*(rv/(rgu+mo/V))*exp(-b*y1*T0/(T0-a*y1))+(rv/(rgu+mo/V))*p2*exp(-b*y1*T0/(T0-a*y1))*y2**2]
         else:
                  return [y2,-g+g*(rv/(rg+mo/V))*exp(-b*y1*T0/(T0-a*y1))-(rv/(rg+mo/V))*p2*exp(-b*y1*T0/(T0-a*y1))*y2**2]
t =arange(0,tz+555,0.1)
y0 = [0.0,0.0]
[y1,y2]=odeint(fun, y0,t, full_output=False).T
plt.title("Подъём, зависание, спуск ЛАЛВ n с жёсткой оболочкой сферической формы  n Объём: %s м3. Масса : %s кг. Подъёмная сила: %s kН. "%(round(V,0),mo,round(0.001*g*rv*V,0)))
plt.plot(t,y1,label='Максимальная высота подъёма: %s км. n Максимальная скорость: % s м/с .n Время зависания %s с.'%(round(max(y1)/1000,2), round(max(y2),2),tz-2*555))
plt.ylabel('Высота в м')
plt.xlabel(' Время в сек.')
plt.legend(loc='best')
plt.grid(True)
plt.show()

Получим:

Динамика вертикального полёта летательного аппарата легче воздуха - 23

Как следует из приведенного графика и листинга программы, для проведения вычислительного эксперимента достаточно ввести необходимые исходные данные.

Математическая модель полёта ЛАЛВ с оболочкой, объём которой изменяется с изменением высоты

К подобным ЛАЛВ относятся стратостаты. Стратостат нельзя полностью надуть гелием, придав ему максимальную подъёмную силу, которая превратит форму его оболочки в шар. Такой шар на большой высоте может лопнуть под действием возросшей разности внутреннего и наружного давлений.

По указанным причинам для расчётов максимально достижимой высоты подъёма вводят два значения его объёма: минимальный Vmin и максимальный Vmax соответственно. С учётом введенных переменных и зависимости плотности воздуха от высоты соотношения для выталкивающей силы Fa и силы тяжести Fт примут вид:

Динамика вертикального полёта летательного аппарата легче воздуха - 24, (4)

Динамика вертикального полёта летательного аппарата легче воздуха - 25, (5)

где: M — масса оболочки и оборудования стратостата;Динамика вертикального полёта летательного аппарата легче воздуха - 26 — плотность гелия.

Приравнивая соотношения (4) и (5), предполагая, что объем оболочки V является функцией от высоты подъёма ЛАЛВ, получим соотношение:

Динамика вертикального полёта летательного аппарата легче воздуха - 27. (6)

Численные значения параметров входящих в соотношение (6) приводятся в листинге для построения графика, который приводится только с указанной целью.

Листинг графика с данными

# -*- coding: utf8 -*-
from numpy import*
from scipy.integrate import odeint
import matplotlib.pyplot as plt
g=9.81# ускорение свободного падения на земле в м/с2.
rv=1.29# плотность атмосферного воздуха в кг/м3.
rg=0.17# плотность гелия в кг/м3.
Vmin=400# начальный объём шара в/м3.
b=0.000128# константа, связанная с плотностью воздуха в 1/м.
c=0.8#коэффициент лобового сопротивления
mo=40#масса в кг
rs=rg+mo/Vmin# суммарная плотность материала стратостата, массы гелия и нагрузки
p1=rv/rs# введенный параметр
h=[(10**-3)*log((rv*w)/(mo+rg*Vmin))*b**-1 for w in arange(1*10**3,1.8*10**5,1000)]
v=[(10**-3)*w for w in arange(1*10**3,1.8*10**5,1000)]
plt.title("Теоретическая зависимость высоты подъёма стратостата n от его максимального объёма")
plt.plot(v,h)
plt.xlabel('Максимальный объём стратостата в тыс. м3')
plt.ylabel(' Максимальная высота подъёма в км.')
plt.grid(True)
plt.show()

Получим:

Динамика вертикального полёта летательного аппарата легче воздуха - 28

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

Выводы:

  1. Получены математические модели двух типов летательных аппаратов легче воздуха, которые позволяют проводить вычислительные экспериенты для оценочного определения параметров таких аппаратов в идеализированных условиях воздушной среды.
  2. Предложенная многоступенчатая схема численного решения системы дифференциальных уравнений позволяет получить вертикальную траекторию летательных аппаратов легче воздуха на этапах подъёма зависания и спуска.

Ссылки

  1. Пара слов про дирижабли
  2. На пути в космос. Стратостаты
  3. Рыжиков Ю.И. Современный Фортран. — СПб.: Корона принт, 2004. – 288 с.

Автор: Scorobey

Источник


* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js