Введение


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

Большое количество публикаций посвящено ЛАЛВ, например, только на нашем ресурсе приведены две очень интересные статьи [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
    03.07.2018 20:55

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

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

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

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

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

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


    1. nickolaym
      04.07.2018 00:56

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


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


  1. Scorobey Автор
    03.07.2018 21:28

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


    1. sshikov
      04.07.2018 19:37

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


      1. technic93
        04.07.2018 19:49
        +1

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


        1. sshikov
          04.07.2018 20:09

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

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

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


          1. technic93
            04.07.2018 20:35
            +1

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


            1. sshikov
              04.07.2018 20:44

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


              1. technic93
                04.07.2018 21:03

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


                1. sshikov
                  04.07.2018 21:05

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


  1. EndUser
    03.07.2018 22:36

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


  1. nickolaym
    04.07.2018 01:06

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


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


  1. maisvendoo
    04.07.2018 07:15

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

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

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


  1. nightwolf_du
    04.07.2018 09:51

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