Блог русскоязычной поддержки Wolfram Mathematica

Падая быстрее скорости звука...

Перевод поста Джейсона Мартинеса (Jason Martinez), Research Programmer, Wolfram Research, Inc.

Нажмите, чтобы скачать данный пост в виде ноутбука Mathematica...

Оригинальный пост: Falling Faster than the Speed of Sound
Перевод сделала Ольга Лавренюк, аспирантка СПбГУЭФ (ФИНЭК)

Общее количество использованных в посте встроенных функций или символов: 64

Список имен используемых встроенных функций и символов в порядке их появления в коде:
ReplaceAll (/.) | WolframAlpha | List ({...}) | Rule (->, ->) | ImageSize | Blank (_) | CompoundExpression (;) | Set (=) | Map (/@) | Function (&) | Part ([[…]]) | Slot (#) | Exp | Most | Cases | RuleDelayed (:>, :->) | LineBox | Pattern (:) | Infinity | Interpolation | InterpolationOrder | Plot | PlotRange | AxesLabel | Style | PlotStyle | Thick | Times (*, ×) | Power (^) | SetDelayed (:=) | Plus (+) | NDSolve | Equal (==) | Derivative (') | Method | Sqrt | PlotLegends | Placed | Bottom | Maximize | Less (<) | NIntegrate | GreaterEqual (>=, ≥) | MaxSteps | Subscript | Dot (.) | D (∂) | Manipulate | Evaluate | Max | Abs | ImagePadding | Italic | Small | Appearance | Delimiter | Column | Button | TrackedSymbols | True | ControlPlacement | SaveDefinitions | ContentSize | Automatic

В один замечательный день в октябре 2012 г., Феликс Баумгартнер (Felix Baumgartner) спрыгнул с высоты 39045 метров или 24.26 миль над Землей из капсулы, поднятой 334-х футовым (101.8 м) стратостатом (размер которого в два раза больше высоты колонны Нельсона в Лондоне на Трафальгарской площади и в 2,5 раза больше диаметра дирижабля Гинденбурга). Сервис Wolfram|Alpha  может проинформировать вас о том, что этот прыжок эквивалентен падению с высоты равной высоте поставленных друг на друга 4,4 гор Эверест, или с высоты, составляющей 93% от длины марафонского забега.

На высоте 24,26 миль над Землей атмосфера очень разряженная и холодная: примерно -14 градусов по Фаренгейту (-25.556°C). Температура, в отличие от атмосферного давления, не зависит линейно от высоты в этих слоях атмосферы. Wolfram|Alpha  может показать нам как она увеличивается и уменьшается в зависимости не только от таких факторов, как уменьшение плотности воздуха при увеличении высоты, но и от поглощения ультрафиолетового излучения в озоновом слое.

WolframAlpha["temperature 39 km", {{"EntrainedTemperaturePlot:AtmosphericLayers", 1}, "Content"}]/.(ImageSize->_) -> (ImageSize->500)

На высоте 39 километров линия горизонта видна на расстоянии примерно 439 миль (706.5 км). В этом слое атмосферы, называемым стратосферой, атмосферное давление составляет всего лишь 3.3 миллибара, что равно примерно 0.33% от значения  атмосферного давления на уровне моря. Иными словами, масса воздуха, находящегося на высоте свыше 39 километров составляет всего лишь 0.32851%  от общей массы воздуха на Земле. С учетом этого можно сказать, что 99.67% массы мировой атмосферы располагается ниже этой отметки. Эта информация была важна Феликсу для решения задачи, стоявшей перед ним — преодоления звукового барьера в свободном падении, так как коэффициент трения о воздух непосредственно связан с атмосферным давлением. Чем меньше воздуха вокруг его тела, тем меньше это трение, а значит, тем большую скорость он может развить при падении с такой большой высоты. Безусловно, по той же причине, для прыжка он был вынужден надеть специальный, снабженный кислородом, костюм, который позволял бы ему дышать и сохранял бы тепло.

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

WolframAlpha["speed of sound 39 km", {{"EntrainedSoundPlot:AtmosphericLayers", 1}, "Content"}]/.(ImageSize->_) -> (ImageSize->500)

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

points = {#[[2]], Exp[#[[1]]]} &/@Most[Cases[WolframAlpha["speed of sound 16 km", {{"EntrainedSoundPlot:AtmosphericLayers", 1}, "Content"}], LineBox[points_] :->points, Infinity][[-1]]] ;

Произведя интерполяцию по полученному набору точек, мы можем получить картину близкую к соответствующему участку высоты, на которой был произведен прыжок Феликса:

Можно заметить, что существует интересный минимум у скорости звука на отрезке, соответствующем середине падения Феликса, идеальный для преодоления звукового барьера.

Итак, насколько же легко преодолеть звуковой барьер в свободном падении? Wolfram|Alpha не может с легкостью ответить на этот вопрос, однако, с помощью системы Mathematica мы сможем узнать ответ на него.

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

0.5 A v² ρ(z) C_d

где v — скорость, C_d — коэффициент сопротивления, A —  максимальная площадь поперечного сечения, а ρ(z) — плотность воздуха на высоте z.

Другим очень важным параметром, который влияет на нашего прыгуна, является ускорение свободного падения g. Это ускорение слабо зависит, в рамках прыжка, от высоты. Оно составляет 9.69 м/c^2 в начале прыжка и 9.831 м/c^2  при приземлении. Это неудивительно, т. к. 24.26 мили — весьма небольшое расстояние относительно радиуса Земли (3956.6 мили).

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

WolframAlpha["air density at 39 km", {{"EntrainedDensityPlot:AtmosphericLayers", 1}, "Content"}]/.(ImageSize->_) -> (ImageSize->500)

Полное уравнение движения имеет вид:

(1)   m v dv/dz = 0.5 A v² ρ(z) C_d - m g(z)

где m — масса Феликса и его снаряжения.

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

Давайте перепишем это уравнение с учетом скорости, как функции от высоты:

m v dv/dz = 0.5 A v² ρ(z) C_d - m g(z)

С набором данных, которыми мы располагаем мы не можем решить его аналитически, однако, мы можем получить численное решение. Сначала определим функции ρ(z) и g(z). Как и до этого, мы можем извлечь данные по плотности воздуха из Wolfram|Alpha:

points2 = {#[[2]], Exp[#[[1]]]} &/@Most[Cases[WolframAlpha["air density at 39 km", {{"EntrainedDensityPlot:AtmosphericLayers", 1}, "Content"}], LineBox[points_] :->points, Infinity][[-1]]] ; airdensity = Interpolation[points2, InterpolationOrder->1] ; (* плотность воздуха в кг/м^3*)

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

gravitationalconstant = 6.67*10^-11; (*Н*м^2/кг^2*) massearth = 5.9721986*10^24; (*кг*) radiusearth = 6367.5*10^3;(*м*) gravity[z_] := gravitationalconstant*massearth/(radiusearth + z)^2

Теперь сделаем несколько предположений относительно остальных физических параметров. Масса Феликса неизвестна, особенно вместе с оборудованием, необходимым для выживания в этих слоях атмосферы, но средняя масса человека составляет примерно 60 кг. Учитывая огромную массу снаряжения, мы можем предположить, что его масса будет около 90 кг или 200 фунтов:

m = 90 (*масса в кг с учетом оборудования*) ;

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

Cd = 1.3(*верхнее значение в нормальных (дозвуковых) пределах для человека*) ;

Площадь поперечного сечения — это еще один сложный момент. Если Феликс был постоянно в вертикальном положении, то площадь поперечного сечения была бы равна квадратному футу или около 0.1 кв. м. (беря средние значения ширины плеч и длины тела). Но в соответствии с отчетами о прыжке, он входил в штопор, что, вероятно, увеличило его поперечное сечение, и, таким образом, добавило порядка 7 кв. футов или 0.7 кв. м. Мы возьмем среднее значение.

A = 0.4(*средняя площадь попереченого сечения в м^2*) ;

Прежде, чем мы посмотрим на то, как его скорость изменилась при спуске, мы должны также узнать, чему равна максимальная скорость свободного падения. В данном случае, уравнение (1) примет вид:

0 = 0.5 A v² ρ(z) C_d - m g(z),

т. к. скорость не должна меняться. Решим это уравнение по переменной v:

v = (m g (z)/0.5 A ρ(z) C_d)^(1/2)

Просуммировав сказанное выше, получим дифференциальное уравнение, описывающее спуск Феликса.

Решим его:

sol = NDSolve[ {m v[z] v '[z] == 0.5 * Cd * A * airdensity[z/1000] * v[z]^2 - m * gravity[z], v[39045] == 0.0001}, v, {z, 0, 39045}, Method->"StiffnessSwitching"]

Plot[{Sqrt[m * gravity[z]/(0.5 * Cd * A * airdensity[z/1000])], soundspeed[z/1000], v[z]/.sol}, {z, 0, 39045}, AxesLabel-> (Style[#, 18] &/@{"Высота, м", "Cкорость, м/с"}), PlotLegends->Placed[{"terminal velocity", "speed of sound", "Felix's velocity"}, Bottom], ImageSize->600, PlotStyle->Thick]

Легко видеть, что скорость Феликса (желтая) быстро увеличивается до тех пор, пока он не достигает или не превышает скорости звука (фиолетовая)  на высоте 30000 метров. Он также ненадолго превышает конечную скорость для сводобного падения с заданной высоты (голубая), увеличивая свою нестабильность по мере того, как атмосфера замедляет его.

Мы можем разделить его скорость на скорость звука для каждого значения высоты, чтобы найти его скорость в числах Маха.

Plot[{v[z]/soundspeed[z/1000]/.sol}, {z, 0, 39045}, AxesLabel-> (Style[#, 18] &/@{"Высота, м", "Число Маха"}), PlotStyle->Thick, ImageSize->600]

Мы также можем использовать систему Mathematica, чтобы найти точку максимума.

maximumz = Maximize[{v[z]/soundspeed[z/1000]/.sol[[1]], 0<z<39045}, z]

По сравнению с официальными отчетами о скоростях на высоте порядка 30000 м, полученное нами значение высоты, на которой достигается максимальное значение скорости меньше на пару километров, вероятно из-за неточного расчета коэффициента сопротивления на сверхзвуковых скоростях. Однако, наше максимальное значение скорости весьма близко к официальному результату — 1.24 числа Маха.

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

На графике ниже показано то, как его скорость менялась бы в зависимости от начальной высоты.

Общее время его падения может быть найдено из функции скорости. Т. к. v = dz/dt, мы можем выразить его как dz/v = dt:

∫1/(m g (z)/0.5 A ρ(z) C_d)^(1/2) dv = ∫dt = t

Возьмем интеграл в пределах от высоты его прыжка до момента раскрытия парашюта на высоте 2500 метров.

NIntegrate[1/(v[z]/.sol[[1]]), {z, 2500, 39045}]

В результате получаем 4 мин. 7.5 с., что близко к реальному времени 4 мин. 19 с. Мы потеряли несколько секунд, что вероятно связано с нашими допущениями относительно коэффициента сопротивления. С небольшой корректировкой мы можем попытаться приблизительно учесть изменение коэффициента Cd, зависящего от скорости движущегося тела.

sol2 = NDSolve[{m v[z] v '[z] ==
     If[v[z] >= soundspeed[z/1000], 1.5, 1] 0.5 * Cd * A * airdensity[z/1000] *
       v[z]^2 - m * gravity[z], v[39045] == 0.0001}, v, {z, 0, 39045},
   Method -> "StiffnessSwitching", MaxSteps -> 100000] ;

Plot[{Sqrt[m * gravity[z]/(0.5 * Cd * A * airdensity[z/1000])],
  soundspeed[z/1000], v[z] /. sol2}, {z, 0, 39045}, AxesLabel-> (Style[#, 18] &/@{"Высота, м", "Cкорость, м/с"}), PlotLegends->Placed[{"terminal velocity", "speed of sound", "Felix's velocity"}, Bottom], ImageSize->600, PlotStyle->Thick]

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

NIntegrate[1/(v[z] /. sol2[[1]]), {z, 2500, 39045}]

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

Возвращаясь к нашему первоначальному решению зададимся вопросом: что было бы, если он не смог бы раскрыть парашют?

NIntegrate[1/(v[z]/.sol[[1]]), {z, 0, 2500}]

v[0]/.sol[[1]]

Он пролетел бы еще приблизительно 44 секунды, вероятно побив рекорд самого длительного свободного падения. Конечно, он ударился бы о Землю со скоростью 53 метра в секунду (примерно 191 км/час). Но, это того не стоит.

Хотя то, что Феликс перевернулся в полете могло легко привести к такому трагическому событию, но, слава Богу, не привело, появляется вопрос: как может потеряться устойчивость полета? Для ответа на этот вопрос мы можем посмотреть на уравнения Эйлера. Если мы рассмотрим Феликса как жесткий прямоугольный параллелепипед, мы сможем найти тензор момента инерции для него:

MomentI = ( {{1/12 (b^2 + c^2), 0, 0}, {0, 1/12 (a^2 + c^2), 0}, {0, 0, 1/12 (a^2 + b^2)}} ) ;

Сделав несколько простых допущений о размерах прыгающего, мы зададим следующие параметры (в метрах):

a = 0.3 ; b = 1.7 ; c = 0.45 ;

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

ω = Ω₁x + Ω₂ y + n z

где x, y, и z — единичные векторы.

Зададим начальные условия:

initialconditions = {t->0, ω_1[t] ->Ω_1, ω_2[t] ->Ω_2, ω_3[t] ->n} ;

Уравнения Эйлера (ω_i — угловые скорости тела):

где:

{I1, I2, I3} == MomentI . {1, 1, 1}

Если ω₁ и ω₂  очень малы по сравнению с частной производной ω₃  в момент времени t, то ω₃ будет почти постоянной. Следовательно, если мы исключим ω₂  из первых двух уравнений, то получим:

\!\( \*SubscriptBox[\(\[PartialD]\), \(t\)]\( \(\*SubscriptBox[\(\[Omega]\), \(1\)]\)[t]\)\)== -(((I2-I3)(I1-I3)Subscript[\[Omega], 3][t]^2)/(I1 I2)) Subscript[\[Omega], 1][t]

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

Насколько нестабильному? Для получения количественной характеристики нестабильности служит следующая интерактивная демонстрация, которая ясно показывает, что начальное вращение вокруг второй оси немедленно вызывает намного более значительное вращение вокруг первой оси, чем наоборот.  

Manipulate[nds = NDSolve[{I1 ∂_tω_1[t] + (I3 - I2) ω_2[t] ω_3[t] == 0, I2 ∂_tω_2[t] + (I1 - I3) ω_1[t] ω_3[t] == 0, I3 ∂_tω_3[t] + (I2 - I1) ω_1[t] ω_2[t] == 0, ω_1[0] == ω10, ω_2[0] == ω20, ω_3[0] == ω30}, { ω_1, ω_2, ω_3}, {t, 0, tMax}] ; Plot[Evaluate[{ ω_1[t], ω_2[t], ω_3[t]} /. nds[[1]]], {t, 0, tMax}, PlotRange -> {All, 3 Max[Abs[{ω10, ω20, ω30}]] {-1, 1}}, ImagePadding -> 30, ImageSize->530, PlotStyle->Thick], "моменты инерции", {{I1, 1, Subscript[Style["I", Italic], 1]}, 0, 5, ImageSize -> Small, Appearance -> "Labeled"}, {{I2, 2, Subscript[Style["I", Italic], 2]}, 0, 5, ImageSize -> Small, Appearance -> "Labeled"}, {{I3, 3, Subscript[Style["I", Italic], 3]}, 0, 5, ImageSize -> Small, Appearance -> "Labeled"}, Delimiter, "начальные угловые скорости", {{ω10, 1, Subscript["ω", 1][0] }, -5, 5, ImageSize -> Small, Appearance -> "Labeled"}, {{ω20, 0.05, Subscript["ω", 2][0] }, -5, 5, ImageSize -> Small, Appearance -> "Labeled"}, {{ω30, 0.05, Subscript["ω", 3][0] }, -5, 5, ImageSize -> Small, Appearance -> "Labeled"}, Delimiter, "интервал времени на котором ищется решение", {{tMax, 100, Subscript[Style["I", Italic], "max"]}, 1, 1000, ImageSize -> Small, Appearance -> "Labeled"}, Delimiter, "начальные повороты вокруг оси\nс небольшими колебаниями", Column[{Button["поворот вокруг оси 1", ω10 = 1 ; ω20 = 0.02 ; ω30 = pe], Button["поворот вокруг оси 2", ω10 = pe ; ω20 = 1 ; ω30 = pe], Button["поворот вокруг оси 3", ω10 = pe ; ω20 = pe ; ω30 = 1]}], {{pe, 0.02, "величина колебаний"}, 0, 1, ImageSize -> Small, Appearance -> "Labeled"}, TrackedSymbols :> True, ControlPlacement -> Bottom, SaveDefinitions -> True, ContentSize-> {550, Automatic}]

Из демонстрации видно, что движение неравномерное, со значительным вращением вдоль каждой оси, а значит прибывать на месте Феликса было бы крайне опасно. К счастью, Феликс смог затормозить вращение и вернуть под контроль ситуацию, в результате чего он смог раскрыть парашют. Это было удивительное зрелище, которое войдет в историю.

Блог принадлежит “Русскоязычной поддержке Wolfram Mathematica"©
При любом использовании материалов блога, ссылка на блог обязательна.
Создано с помощью Wolfram Mathematica 9