МЕХАНИКА ЖИДКОСТИ И ГАЗА № 3 • 2013
УДК 533.601.1:517.95
© 2013 г. В. А. АНДРУЩЕНКО, И. В. МУРАШКИН, Ю. Д. ШЕВЕЛЕВ
ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ О ВЗРЫВЕ В АТМОСФЕРАХ ПЛАНЕТ В ПЕРЕМЕННЫХ ЛАГРАНЖА
Рассмотрена задача о взрыве малой мощности в атмосферах Земли и Юпитера в переменных Лагранжа. В ходе вычислительного эксперимента выявлены особенности течения внутри возмущенной области, охваченной ударной волной. Помимо сингулярности типа узел, заложенной в начальных условиях задачи, в этой области возникают и другие особые точки — фокусы, седла, узлы. Определено влияние определяющих параметров на структуру течения.
Ключевые слова: астероидно-кометная опасность, точечный взрыв, ударная волна, переменные Лагранжа, особые точки.
В настоящее время задачи теории точечного взрыва вновь выдвинулись в ряд задач первостепенной значимости. Это связано с актуальной проблемой астероидно-комет-ной опасности — возможного входа космических объектов в атмосферу Земли (и других планет) и их взрыва в ней или на поверхности с катастрофическими последствиями. Так, в течение последних 100 лет наша планета подвергалась ударам четырех небесных тел: Тунгусского в 1908 г., Бразильского в 1930 г., Сухотэ-Алиньского в 1947 г. и Витимского в 2002 г. [1—4], которые по счастливой случайности произошли в практически безлюдной местности. А кроме того, десятки опасных тел пролетели в непосредственной близости от Земли [1, 2, 5]. В большинстве разработанных на настоящее время моделях разрушения метеороидов под взрывом понимался переход их начальной энергии в другие виды энергии: торможение, нагрев, светимость и пр. [6]. Но в некоторых работах были рассмотрены модели мгновенного дробления крупных метеорных тел на множество мелких фрагментов, в которых предполагалось, что за счет более быстрого торможения мелких фрагментов головная ударная волна трансформировалась во взрывную волну, производя эффект воздушного взрыва [7,8]. Этой же концепции — взрыва в атмосфере Земли придерживались авторы статей [9—13] и в атмосфере Юпитера [14].
1. Постановка задачи. Исходные уравнения и краевые условия. Методика расчета. В [14] было предложено моделирование взрыва фрагмента кометы в атмосфере Юпитера посредством цилиндрического взрыва в экспоненциальной (изотермической) атмосфере (на начальном этапе — автомодельным решением о сильном взрыве [15]). Действительно, в диапазоне высот 50 км < к < 320 км атмосфера Юпитера изотермична (фиг. 1, а [16]), а распределения давления и плотности достаточно хорошо аппроксимируется экспоненциальной зависимостью
где р0, р0 — давление и плотность на высоте взрыва, а А — масштаб неоднородности атмосферы (см. ниже).
(1.1)
Фиг. 1. Зависимость давления (мбар) от температуры К в атмосферах Юпитера (а): 1 — тропопауза 50 км, 2 — стратопауза (320 км), 3 — верхняя граница (1000 км); и Земли (б): 1 — тропопауза 11 км, 2 — стратопауза 26 км, верхняя граница 40 км
Фиг. 2. Расчетная область
В настоящей работе использовано решение в приближении сферического взрыва, как моделирующего мгновенное разрушение малого космического объекта с формированием взрывной волны (на начальном этапе — автомодельное решение о сильном точечном взрыве [15]), как это было предложено в [9—11] для взрыва в атмосфере Земли при прочих аналогичных предположениях относительно атмосферы Юпитера в указанном выше диапазоне высот.
За исходную систему дифференциальных уравнений выбрана система уравнений газовой динамики в лагранжевых переменных £ и 0. Координата частицы £ определена моментом времени, когда фронт ударной волны доходит до этой частицы, а за координату 0 принят угол, образованный нормалью к фронту ударной волны в момент прохождения через данную частицу с вертикальным направлением (осью г). В силу такого выбора, ударная волна в лагранжевых координатах £, 0, г задана уравнением £ = t, а область возмущенного газа внутри волны отображена в область t > t0, 0 < £ < t (фиг. 2)
Исходная система уравнений в лагранжевых переменных в безразмерном виде записывалась следующим образом:
-(-Ü z • + £'«) (Р '
-(-is*+t*) с- a
ди
дг
ди
дг
др = -ур(дгд£-дгди+д1ди-д1ди}-1Ри (1.2)
дг I [д^де двд^ двд^ д^дв^ г
д (р1г) = о, дГ = и, = и
дг дг дг
I = г^г 0 - г0 г ^
Искомые функции: г, z — эйлеровы координаты частицы, и, и — составляющие скорости, р, р — плотность и давление возмущенного газа. Здесь у — показатель адиабаты, = £р0Л4 /Е0 — определяющий безразмерный параметр гравитации, g — ускорение свободного падения, Е0 — энерговыделение при взрыве.
Граничные условия следующие. На фронте ударной волны Г! (фиг. 2) при £ = г это условия Ренкина—Гюгонио
:— N cos a \1-у |, и =— N sin a íl-y APr Y +1 l N2) Y +1 l N2
Pi 11 + ^7I expI-—I (1.3)
:X±I|l + 2L. Ai exp |-— Y-11 Y-1 N ) \ A,
Y -1 ap —
P1 =N I 1 |exp|-
Y +1 l 2 N ) V A
Здесь N — скорость ударной волны, a — угол между нормалью к фронту и осью r,
Ap = (po /E0) A3 — определяющий безразмерный параметр противодавления. Величины с индексом 1 — значения газодинамических величин на фронте ударной волны.
Эйлеровы координаты точек на ударной волне определены геометрическими соотношениями
dr лт • n , dN n d— лт n dN • Л
— = N sin 0+--cos 0, — = N cos 0--sin 0
dt 50 dt d0
В центре взрыва £ = 0 выполнены условия симметрии
u = 0, r = 0, d-P = йР = ди = o
д r д r д r
Расчетная область разбивалась на две подобласти: малую внутреннюю G0 и протяженную внешнюю G1 (фиг. 2). Для удобства разностного счета подобласть G1 преобразовывалась в прямоугольник посредством перехода к новой переменной
Ь = d t "У (1.4)
При этом внешняя граница внутренней области ^(t) преобразовывалась в координатную линию А = 1. Лучи, покрывавшие расчетную область, могли смещаться со временем, что позволяло обеспечить "более правильную" конфигурацию ячейки сетки в физическом пространстве. Кроме того, над лучевой координатой £ осуществлено преобразование сжатия, которое ставило в соответствие равномерной сетке в расчетных координатах существенно неравномерную сетку в физической плоскости (сгустив уз-
лы сетки вблизи ударной волны, где искомые газодинамические величины имели большие градиенты), не изменив прямоугольного вида подобласти Преобразование сжатия зависило от управляющего параметра р в (1.4), изменяя который можно было подбирать желаемое распределение узлов сетки вдоль луча и вносить в нее коррективы в ходе расчета.
В малой центральной подобласти О0 (фиг. 2) для нахождения искомых функций использовались следующие асимптотические формулы:
Нижние индексы 0—2 отвечают точкам с лагранжевыми координатами X = 0, 1/и. 2/п, где п — количество расчетных точек по оси X.
Во внешней протяженной подобласти ^ (фиг. 2) искомые функции рассчитывались с помощью явного двухшагового разностного метода типа предиктор—корректор [17].
2. Тестирование методики. Результаты. Для тестирования сконструированного алгоритма и написанной программы просчитывалась задача о сильном взрыве в отсутствие противодавления и гравитации (Ар = А& = 0) в однородной атмосфере (А ^ да), имеющая аналитическое решение [15]. Результаты численного расчета совпадали с точным решением с очень высокой степенью точности. Результаты расчетов вариантов с противодавлением и гравитацией, проведенные на основной разностной сетке с числом узлов 1600 (40 х 40) и двойной — 6400 (80 х 80), совпали до многих знаков после запятой. Контроль основных расчетов далее велся посредством слежения за выполнением законов сохранения массы и энергии. Так, полная энергия всей возмущенной массы газа, вычисленная по рассчитанным значениям функций в любой момент времени ?, должна равняться сумме энергии взрыва и полной энергии невозмущенного газа в объеме, охваченном ударной волной в текущий момент времени ?
Здесь для проведенных ниже расчетов ю(у) = 0.851 при у = 1.4.
Поставленная выше задача о взрыве зависит от трех безразмерных параметров, у, Ар, А& и одного размерного параметра А, отвечающего за скорость экспоненциального изменения плотности и давления в изотермической атмосфере (1.1).
Предыдущие расчеты [13] показали, что решение для газа (1 < у < 2) качественно слабо зависит от у, при меньших значениях параметра газ только в большей степени концентрируется у ударного фронта, что значительно увеличивает прифронтовые градиенты искомых функций. Кроме того, было установлено, что в отсутствие хотя бы одного физического фактора — гравитации (Ая = 0) или противодавления (Ар = 0) в ходе решения в поле мгновенных линий тока имела место только одна особая точка типа узел — сингулярность, внесенная еще в начальные условия выбором автомодельного решения Л.И. Седова. В ходе эволюции течения внутри возмущенной области при таких условиях происходит только перемещение этой особой точки вдоль линии
где Ь = (1/2)
2(у-1)/5(уР)
(2.1)
симметрии — оси %. Здесь следует отметить, что безразмерный параметр Ая присутствует только в одном из исходных уравнений (второе уравнение в (1.2)), а параметр Ар — только в граничных условиях (1.3). И лишь при Аё ^ 0 и Ар ^ 0 эволюция течения сопровождается формированием дополнительных особых точек типа фокусов, седел и узлов [13].
В настоящей работе на первом этапе проводился численный эксперимент для выяснения роли единственного размерного параметра — А на характер эволюции течения в возмущенной области при фиксированных численных значениях всех безразмерных параметров у, А& и Ар. Для этого рассматриваются взрывы в атмосферах Земли и Юпитера, где ДЕ = 6.4 км и А 7 = 23.5 км. Для равенства параметров Ар = 5, Ая = 5.2 при у = 1.4 выбираются взрывы с начальным энерговыделением Е0/ = 1 Мт на высоте И01 = 90 км на Юпитере и Е 0Е = 82 Кт на высоте Н0Е = 18 км на Земле.
Опишем как проистекает эволюция структуры течения внутри возмущенной области в этих случаях в атмосферах Земли и Юпитера.
Высота взрыва в атмосфере Земли выбирается равной 18 км для того, чтобы она попала в интервал изотермичности ее реальной атмосферы — 11 км < к < 26 км (фиг. 1, б) [18]. Расчет начинается с радиуса взрыва Д = 200 м (при этом скорость газа на фронте ударной волны равна 6680 м/с, давление 6.75 МПа). Усл
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.