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

Введение


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

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

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

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

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

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


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





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



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

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

, (1)

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

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

, (2)

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

;
;

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

, (3)

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


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



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



где: 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()


Получим:



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

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


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

Подъём;
Зависание ;
Приземление.

Листинг программы
# -*- 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()


Получим:



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

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


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

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

, (4)

, (5)

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

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

. (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()


Получим:



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

Выводы:

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

Ссылки


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

Вы можете помочь и перевести немного средств на развитие сайта



Комментарии (14):

  1. sshikov
    /#18843147

    >по закону Миндалева Клапейрона

    вот Дмитрий Иванович на вас обиделся :)

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

    Это вы какое сопротивление хотите учесть, то которое при подъеме со скоростями порядка метров в секунду? Гм. Есть подозрение, что им спокойно можно пренебречь.

    >площадь Миделя

    Вообще-то мидель это не фамилия.

    • nickolaym
      /#18844007

      товарищ Мидель Кастро :)))


      Менделеевый Клайперон продолжает недоумевать без тире.

  2. Scorobey
    /#18843257

    За указания на опечатки спасибо а вот с арифметикой хуже. Сила сопротивления 0.5*с*s*v2=0.5*0.4*3.14*625*40=15,7 kН. Можете не опасаться. Спасибо за комментарий.

    • sshikov
      /#18848593

      Гм. Плотность воздуха 1,2754 кг/м?. Что-то я не вижу у вас такого числа, даже близко. Вижу, что шарик видимо, судя по Cx, размером 25 метров, вижу что скорость где-то 6-7 метров в секунду (40 это наверное квадрат?). Осталось понять, каким образом этот шарик радиусом 25 метров? набрал 6 метров в секунду? Обычно такие стратостаты так быстро не летают.

      • technic93
        /#18848625 / +1

        Сила трения как раз и ограничивает максимальную скорость.

        • sshikov
          /#18848691

          Во-первых, не трения. Сопротивление трения для шарика — всего-лишь примерно 10%.

          Еще раз поясню — мне лично кажется, что получить такую скорость, чтобы аэродинамическое сопротивление ее ограничивало, за счет подъемной силы, невозможно практически. Ну то есть, понятно что если ваш воздушный шарик попал в ураган, и сила ветра скажем метров 30 в секунду — то да, сопротивление воздуха нифига игнорировать не получится. Для падения тоже самое — там ускорение 9.81 делает свое дело.

          Подъемная сила — она же тоже не бесконечна, максимум что мы можем — наполнить шарик водородом или гелием. Ну и практически я исхожу из того, что гелиевые воздушные шарики поднимаются на 100 метров вовсе не за 10 секунд, а скажем за минуту — т.е. скорости подъема даже 10 метров в секунду — это на глаз какая-то фантастика.

          • technic93
            /#18848767 / +1

            Т.е. по вашему шарик будет все время ускоряться?

            • sshikov
              /#18848797

              С чего вы взяли? Я же ясно написал, что будет равновесная скорость, потому что тут есть единственная сила, которая пропорциональна ее квадрату (а две другие более-менее постоянны). И я именно хочу понять, какова оценка этой скорости, и из чего она вытекает.

              • technic93
                /#18848843

                В вашем первом комментарии вы предлагаете ей принебречь.

                • sshikov
                  /#18848851

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

  3. EndUser
    /#18843533

    Опять без итерационного метода?
    «Домашка»?

  4. nickolaym
    /#18844027

    Как внезапно, листинг на питоне, библиография — "Современный Фортран" :)


    Я так понимаю, из книги по фортрану были взяты формулы метеорологии и аэродинамики?

  5. maisvendoo
    /#18844401

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

    Каких, к черту, линейных, когда одна из фазовых координат (скорость) в квадрате, а другая — (высота) входит как аргумент функции exp(x).

    Не путайте линейность с формой Коши. И вообще, матчасть получше изучите…

  6. nightwolf_du
    /#18844865

    Возможно, глупый вопрос из-за незнания мат.части — но что мешает для шариков использовать ту же комбинацию из акселерометров-гироскопов-магнетометра-барометра что и для дронов?
    habr.com/post/137595