Безразмерный воздушный шар. Утилитарная магия анализа размерностей05.07.2018 08:03

ezw2yecjcancuztlb2fzelxeleu.gif

На написание этой небольшой заметки меня натолкнула недавно опубликованная на Хабре статья Динамика вертикального полёта летательного аппарата легче воздуха. Захотелось написать комментарий, но он быстро перерос во что-то большее и, как кажется, более полезное.

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

Можно сделать так, чтобы модель и расчёты стали универсально полезными не для какого-то конкретного шара, а для широкого круга задач. Можно обеспечить оптимальную точность вычислений при численном интегрировании дифференциального уравнения. Можно избавиться от необходимости вручную задавать пределы интегрирования и шаг при расчёте в широком диапазоне параметров. Наконец, можно многое рассказать о динамике полёта нашего шара и без численного решения. И для всего этого служит один давний приём, верный и надёжный, когда-то обязательный при любых расчётах на ЭВМ и до их появления, а сейчас факультативный и часто относимый к магии и искусству — приведение уравнений к безразмерному виду и собственным масштабам. Воспользуюсь задачей о воздухоплавании, как примером и покажу, насколько более осмысленным и изящным становится анализ задачи, при использовании этой техники. А потом объясню почему это может быть важным для программистов, и отчего эта статья попала в хаб «Функциональное программирование».

Суть приведения уравнений к собственным маштабам состоит в отыскании и введении таких единиц измерения времени, расстояний, масс, токов и прочих размерных переменных, в которых задача принимает наиболее простой и элегантный вид. Такие единицы называются характерными или собственными масштабами задачи. При этом существенно сокращается число параметров, влияющих на её решение, а если повезёт наткнуться на автомодельное решение, то параметры могут исчезнуть вовсе!

Обычно, с этим приёмом знакомятся, когда изучают динамику гармонического осциллятора, приводя уравнение движения к каноническому виду. Действительно, каким бы ни был осциллятор (большим или маленьким, низко- или высокочастотным, механическим или электромагнитным), если измерять время не в секундах, а в периодах его свободных колебаний, а амплитуду не в метрах, а, скажем, в величине начального отклонения, то уравнение для любого линейного осциллятора будет иметь один и тот же вид:

$x'' + 2\zeta x' + x = 0$

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

Также дело обстоит и при анализе кривых второго порядка. Рассмотрев параболу, гиперболу или эллипс в их канонической (самой компактной) форме, мы узнали все их замечательные свойства. После чего, эти знания распространяются на любую нетривиальную квадрику. Более того, этот же приём используется при анализе уравнений матфизики второго порядка (в частных производных).

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

  1. Вводим формальные масштабные множители для переменных;
  2. делим все уравнение на какую-либо величину, имеющую размерность его слагаемых;
  3. из получившихся безразмерных комплексов, включающих масштабные множители, выбираем столько, сколько множителей у нас есть, и приравниваем их единице;
  4. решаем получившиеся уравнения, выражая масштабные множители через параметры задачи, а оставшиеся комплексы используем, как критерии подобия.


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

Давайте обратимся к исходному уравнению движения, приведённому в статье о летательных аппаратах, и покажем на его примере, как можно выполнить его преобразование к безразмерному виду в собственных масштабах. Уравнение имеет вид:

$ m \frac{d^2h}{dt^2} = -m g + g W \rho_0 e^{-b h} - \frac{1}{2}c S \rho_0 e^{-b h} \operatorname{sign}\left(\frac{dh}{dt}\right)\left(\frac{dh}{dt}\right)^2, $

с начальными условиями

$h(0) = h'(0) = 0.$

Здесь $h$ — высота подъема шара, $m$ — масса всего летательного аппарата с грузом, $g$ — ускорение свободного падения, $W$ — объем газа в шаре, $c$ — коэффициент лобового сопротивления, $S$ — характерная площадь сопротивления, $ρ_0$ — плотность воздуха на нулевой высоте, $b$ — коэффициент в распределении Больцмана.

Первый шаг. Введём формальные масштабы для времени и расстояний:

$ h = h_0 y,\quad t = t_0 \tau \qquad (1) $

и перепишем уравнение движения, используя штрихи, для обозначения производных:

$ m \frac{h_0}{t_0^2} y'' = -m g + g W \rho_0 e^{-b h_0 y}\left[1 - \frac{c S}{2W g}\frac{ h_0^2}{t_0^2} \operatorname{sign}(y')(y')^2\right]. $

Начальные условия в нашем случае тривиальны, поэтому их можно опустить, но вообще, надо переписывать и их.

Второй шаг, собственно, обезразмеривание уравнения. Все слагаемые в нём имеют размерность силы и мы вольны поделить их на любую силу. Обычно делят на силу инерции — на множитель при второй производной от расстояния. В таком случае, в качестве параметров задачи, обычно получаются хорошо знакомые критерии динамического подобия, такие как число Рейнольдса или Эйлера. Но в нашей задаче мы все силы, входящие в задачу, отнесём к силе тяжести $mg$, и вот почему. Нас интересует, главным образом, положение статического равновесия — предельная высота полёта шара при заданной грузоподьъёмности и время перехода к нему. А статическое равновесие не зависит от инерционных свойств системы, зато напрямую зависит от гравитации. Итак, делим и сокращаем, что возможно:

$\frac{h_0}{gt_0^2} y'' = -1 + \frac{W \rho_0}{m} e^{-b h_0 y}\left[1 - \frac{c S}{2W g}\frac{ h_0^2}{t_0^2} \operatorname{sign}(y')(y')^2\right].$

Всё, уравнение обезразмерено, теперь все переменные и все слагаемые в нём — это просто числа. И теперь мы свободны выбирать такие масштабы длины и времени, чтобы максимально сократить число параметров задачи. У нас есть две неизвестные $h_0$ и $t_0$ и четыре безразмерных комплекса:

$ \frac{h_0}{gt_0^2},\quad \frac{W \rho_0}{m},\quad b h_0,\quad \frac{c S}{2W g}\frac{ h_0^2}{t_0^2} $

значит, мы можем заставить два из них исчезнуть, просто приравняв их единице. Причём второй комплекс не зависит от масштабных множителей, и нам надо выбрать для уничтожения два оставшихся. И опять мы обратимся к аргументу об интересующем нас статическом равновесии. Положение равновесия определется равенством нулю ускорения, то есть, нам будет неважно какой множитель окажется при второй производной от высоты и можно заняться минимизацией числа параметров в правой части уравнения. Таким образом, мы приходим к следующему решению: масштабные множители выберем так, чтобы третий и четвёртый из перечисленных нами комплексов исчезли из уравнения. Для этого полагаем

$ b h_0 = 1,\quad \frac{c S}{2W g}\frac{ h_0^2}{t_0^2} = 1 $

и получаем собственные масштабы задачи:

$ h_0 = \frac{1}{b},\quad t_0 = \frac{1}{b}\sqrt{\frac{c S}{2W g}}.\qquad (2) $

Оставшиеся параметры обозначим следующим образом:

$ \frac{h_0}{gt_0^2} = \frac{1}{\gamma},\quad \frac{W \rho_0}{m} = B. $

Это, наверное, самый магический этап наших преобразований, который ближе всего к искусству. Но опыт и некоторые рассуждения помогают в выборе правильных масштабов. А потом этот же опыт позволяет видеть подобные явления в разных областях науки и угадывать характерные свойства сложных систем.

Я позволю себе немного отложить анализ физического смысла полученных нами масштабов и параметров, хоть это тоже очень приятно и интересно, но сразу перейду к самому вкусному — к безразмерному уравнению в собственных масштабах:

$ \frac{1}{\gamma} y'' = B e^{-y}\left(1 - \operatorname{sign}(y')(y')^2\right)-1,\quad y(0) = y'(0) = 0.\qquad (3)$

Вот! Такое уравнение я готов решать, анализировать и численно интегрировать. А когда решу — легко перейду от собственных масштабов к метрам и секундам, закодированным в соотношениях (1) и (2).

Итак, начнём. Найдём, для начала, где наш обобщённый воздушный шар остановится. Для этого приравняем нулю как скорость, так и ускорение в уравнении (3):

$ 0 = B e^{-y^*} - 1, $

откуда получим

$ y^* = \ln B. $

Сразу делаем вывод: для того, чтобы предельная высота существовала, требуется выполнение условия $B>1$» />. Далее, оценим, за какое время шар достигнет этой высоты. Для этого приравняем нулю ускорение и решим задачу Коши первого порядка: </p> <p><img src=

Оно решается аналитически методом разделения переменных, и сводится к интегралу:

$ \tau^* = \int_0^{y^*} \frac{B dy}{\sqrt{B-e^y}} = 2 \operatorname{arth}\left(\sqrt{1 - \frac{1}{B}}\right). $

Смотрите, вот они — плоды правильно выбранных масштабов! И предельная высота, и характерное время её достижения выражаются только лишь через параметр $B$. Это значит, что мы можем зафиксировав какое-либо значение $B>1$» />, построить ряд графиков решения уравнения (3) для разных значений <img src= и тем самым сразу описать все варианты решения задачи для любых значений$B$!

Вот это универсальное однопараметрическое семейство решений, полученное численным интегрированием задачи (3) для различных значений $\gamma$:
soacp1f33kwxgw6ggeuy3hzcapw.png

После достижения равновесной высоты, воздушный шар какое-то время совершает затухающие колебания, однако характерное время достижения этой высоты оценено корректно. От критерия $\gamma$ зависят лишь период колебаний и время их затухания. Чем этот параметр больше, тем «жёстче» наша система.

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

Теперь посмотрим, что же такое у нас получилось в виде параметров и переменных. Заменим массу шара $m$ на сумму массы газа $W\rho_g$ и массы полезного груза $M$:

$m = W \rho_g + M = W \rho_g (1 + \alpha),$

где $\alpha = \frac{M}{W\rho_g}$ — относительная грузоподьёмность летательного аппарата. Кроме того, примем шар сферическим и выразим его площадь и объём через радиус $R$. В таком представлении мы получаем параметры задачи:

$B = \frac{\rho_0}{\rho_g (1 + \alpha)},\\ \gamma = \frac{3 c B}{8 R b}.$

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

Масштабы выражаются таким образом:

$h_0 = 1/b\\ \tau_0 = \frac{1}{2 b}\sqrt{\frac{3 c B}{2 g R}} = \sqrt{\frac{\gamma}{g b}}$

Масштаб расстояния определяется только градиентом плотности. Что совершенно справедливо, ведь именно из-за этого градиента подъём шара, вообще, где-либо останавливается. Характерное время включило в себя динамические величины — ускорение свободного падения, массовые соотношения и сопротивление воздуха.

Характерные высота и время подлёта шара

$h^* = y^* h_0 = \frac{1}{b}\ln B,\\ t^* = \tau^* t_0 = 2\sqrt{\frac{\gamma}{g b}} \operatorname{arth}\left(\sqrt{1 - \frac{1}{B}}\right)$

также выражаются через основные параметры и масштабы задачи.

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

$B = 4.57\\ \gamma = 686\\ h_0 = 8000\ м\\ \tau_0 = 12.47\ мин\\ h^* = 12166.6\ м\\ t^* = 13.72\ мин$

Подведём итог. Немного поработав на бумаге перед тем, как скармливать задачу машине, мы смогли получить более внятный, универсальный результат и подметили самые важные свойства решения. При этом, мы не потеряли «модульности» решения. Вот что я имею в виду. Мы сформулировали задачу в самых общих чертах, а потом, уже получив решение, начали добавлять новые «фичи», например, разделили понятие полезного груза и массы поднимающего его газа. Это усложнение повлияло на значения масштабов, но не изменило характера решения.

Если впоследствии надо будет добавить учёт температурного градиента, выражение в экспоненте усложнится, но его суть и, главное, масштаб не изменятся:

$\exp\left(\frac{-b h T_0}{T_0 - a h}\right)$

после обезразмеривания будет выглядеть так:

$\exp\left(\frac{-y}{1 - \psi y}\right), $

где новый параметр $\psi = \frac{a}{bT_0}$ показывает насколько соотносятся между собой градиенты плотности и температуры. Так как это параметры атмосферы, от нас не зависящие, $\psi$ — это константа, равная $0.17$. Величина константы, кстати, показывает роль эффекта температурного градиента, он не велик, но и не пренебрежимо мал. Семейство графиков немного изменится, но, главное, останется однопараметрическим.

У использования собственных масштабов задачи при расчётах есть ещё один существенный плюс: переменные при этом, обычно, принимают «умеренные» значения, то есть, близкие к единице. Это очень полезно для вычислений с плавающей точкой: не теряется точность при операциях над значениями, сильно различающимися по порядку величины. Кроме того, появляется возможность корректно сравнивать роли отдельных эффектов в задаче, сравнивая величины соответствующих им критериев подобия или значений переменных. Например, значение $B=4.5$ говорит о том, во сколько раз силы взаимодействия с атмосферой больше силы тяжести. А квадрат характерной безразмерной скорости подьёма $(y^*/\tau^*)^2 \sim 0.3$ показывает степень значимости сопротивления воздуха при вертикальном движении летательного аппарата.

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

А что было бы, если бы мы выбрали другие комплексы, или отнесли все силы не к силе тяжести, а, скажем, к инерции? В уравнении осталось бы два параметра, но семейство кривых стало бы двухпараметрическим, и его нельзя было бы показать на одной диаграмме. Признаюсь, сначала у меня именно так и вышло, в качестве масштаба длины я получил радиус шара. Но поиграв с этими кривыми, я увидел их геометрическое подобие (так во всех параболах видишь одну), а дальнейший беглый анализ получившихся критериев и уравнений сам натолкнул меня на скрытую симметрию и подсказал каким должен быть «канонический вид» наших уравнений. Это, правда, красивое занятие! Задача сама начинает о себе рассказывать. Подобное удовольствие я испытываю выстраивая систему типов и структур данных в программе на Хаскеле или на C#: когда архитектура программы соответствует внутренней структуре задачи, всё становится удивительно естественным, элегантным, «автомагически» отрабатываются особые случаи и уменьшается число слоёв абстракции.

Размерности физических величин сами по себе очень интересны. Они образуют линейное пространство, и поиск критериев подобия можно свести к задаче поиска ядра в пространстве размерностей, тем самым формализуя этот процесс. Они, в какой-то мере, играют роль типов в физических вычислениях. Как компилятор использует статическую типизацию для проверки корректности программы, так и физик использует размерности для верификации своих выкладок и результатов. Как в строго типизированном чистом функциональном языке (например, в Хаскеле) можно вывести код функции из её типа, таким же образом, исходя из размерности физических величин, можно конструировать безразмерные комплексы, характерные величины в физических системах и получать чрезвычайно полезные и универсальные результаты. Тому есть множество примеров в механике, газо- и термодинамике, квантовой механике и т.п. Я рекомендую познакомиться с замечательной работой, в которой приводится куча красивых примеров применения анализа размерностей к задачам от теоремы Пифагора, до колебаний звёзд и Релеевского рассеяния в небе. В этой работе цитируются слова Джона Уилера — учителя Ричарда Фейнмана, которые получили имя «правила Уилера»:

«Никогда не начинай вычислений, пока не знаешь ответа. Каждому вычислению предпосылай оценочный расчет: привлеки простые физические соображения (симметрию! инвариантность! сохранение!) до того, как начинать подробный вывод; продумай возможные ответы на каждую загадку. Будь смелее: ведь никому нет дела до того, что именно ты предположил. Поэтому делай предположения быстро, интуитивно. Удачные предположения укрепляют эту интуицию. Ошибочные предположения дают полезную встряску».


Этот совет для программистов звучит очень знакомо: не начинайте писать код функции, не определившись с её сигнатурой (типом) и её поведением (тестами); прежде написания кода следует продумать типы, структуры данных и их взаимоотношения. Именно этим сильны ООП и функциональное программирование — там эти принципы не просто работают, они, при искусном использовании, сами подсказывают наиболее естественные и изящные решения, а в случае с ФП, делают это на глубоком математическом уровне, давая возможность доказывать свойства программ и поручать часть работы по выводу свойств программы компилятору.

Изучайте математику, программируйте красиво и получайте удовольствие!

© Habrahabr.ru