CC BY
50
Вычислительные технологии
Том 13, Специальный выпуск 2, 2008
Моделирование падения в атмосфере остатков ракетного топлива*
Ю.Н. Мороков
Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: quant@ict.nsc.ru
A further development of the model of a falling ensemble of evaporating rocket fuel drops in the atmosphere is presented. The new model allows considering an emission of fuel at any moment, starting from the moment of a rocket stage separation at high altitudes. Consecutive dynamic modelling of trajectories of falling fuel drops in the atmosphere is presented. Results of calculations illustrate that an emission of the fuel residues from a rocket stage at high attitudes when the stage has still sufficiently high speed, solves the problem of elimination of the ecological risks related to pollution of the terrestrial surface by the rocket fuel residues. Numerical modelling of time evolution of a vapour cloud, formed due to the evaporation of falling drops of fuel in the atmosphere is also considered. The problem has its own specifics as both the localized moving source (a falling rocket stage) and the distributed moving source (falling evaporating drops) are taken into account. An opportunity of the application of a universal method of construction of adaptive grids, based on the numerical solution of the inverted Beltrami and diffusion equations is considered in this problem.
Отработавшие вторые ступени ракет-носителей "Протон-М" отделяются на высотах порядка 130 км и разрушаются при входе в плотные слои атмосферы па высотах 30...40 км. Остатки ракетного топлива (несимметричного диметилгидразина — НДМГ) попадают в атмосферу после разгерметизации падающих ступеней. Это токсичное вещество проявляет канцерогенные и мутагенные свойства и относится к суперэкотокси-кантам. НДМГ растворим в воде, легко испаряется и может загрязнять воздух, почву, наземные и подземные воды.
Остаточное количество ракетного топлива в баке второй ступени ракеты-носителя "Протон-М" составляет около 300 кг. Выброс топлива из падающей ступени может происходить вдоль всей траектории ее падения. При торможении в атмосфере топливо будет дробиться на капли с диаметром до 6 мм, которые при падении будут испаряться.
Ранее нами проведены расчеты загрязнения поверхности земли ракетным топливом при падении отделяющихся частей ракет-носителей "Протон" для девяноста конкретных эпизодов 11 3| в трех районах падения на территории Республики Алтай и Восточного Казахстана. В основу использованной в [1-3] физической модели падения в воздухе ансамбля испаряющихся капель НДМГ положена модель, разработанная ранее Э.Л. Александровым [4|.
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 06-01-08009, № 07-01-00336 и № 08-01-00264).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.
В настоящей работе представлено дальнейшее развитие модели. В предыдущих расчетах рассматривалось падение в атмосфере капель ракетного топлива с высот не более 30 км. При этом полагалось, что капли падают с квазиравновесными установившимися скоростями. В новой модели эти ограничения снимаются и появляется возможность рассмотрения целенаправленного раннего выброса топлива па больших высотах с целью минимизации экологических рисков на земной поверхности и в нижних слоях атмосферы. Модель включает последовательное динамическое моделирование траекторий падения испаряющихся капель топлива в атмосфере. Система динамических уравнений, описывающих движение испаряющейся капли массы т(Ь), в сферической системе координат в инерциальной системе отсчета имеет следующий вид:
dr d0 dу
dt r dt ' dt
uu u,^ ,
Vr,— =We,— =Wlp, (1)
^ = 1 (/, + U) + rwl + r cos 2(d)wl, (2)
- f'ce - cos(û) sin[9)w2 - -vrwe, (3)
dt mr v r
fcv-2wv(^-wetg(e)), (4)
dwv 1 fVr
dt mr cos(e) V r
(5)
dm
— = 2тг rkp, dt
ark\l + 2(1-a)D
Здесь r — расстояние до центра Земли; в и у — широта и долгота (в невращающейся системе отсчета) ; vr — радиальная скорость капли; wq и wv — соответствующие угловые скорости; fg — ньютоновская сила притяжения капли к Земле; fcr, fcQ и fcip — проекции на соответствующие координатные орты силы сопротивления, действующей со стороны воздуха на каплю; rk — радиус капли (предполагается, что она имеет сферическую форму) ; pn — плотность насыщающих паров НДМГ при текущей температуре поверхности капли Ts (К) ; R — универсальная газовая постоянная; ^ — молярная масса НДМГ; D — коэффициент диффузии молекул НДМГ в воздухе [1].
Параметр 0 < а < 1 в уравнении (5) используется для линейного интерполирования между двумя предельными механизмами отвода молекул НДМГ от капли. Первое слагаемое в скобке учитывает только вылет молекул с поверхности капли, но не их обратную конденсацию на капле. Оно работает для очень разреженной атмосферы при больших числах Кнудсена. Второе слагаемое становится существенным, когда скорость изменения массы капли лимитируется диффузионным отводом молекул НДМГ в атмосферу [1, 4]. В предварительных расчетах мы использовали это простое линейное интерполирование для переходной области высот 30 < z < 80 км.
В рассматриваемой модели мы пока не учитываем фазовый переход НДМГ из жидкой в твердую фазу (НДМГ кристаллизуется при температуре -58 °С). В обсуждаемом ниже примере численного моделирования падения капли, образовавшейся сразу же по-
Ts
диапазоне от -115.5 до -58.3 °С, что соответствует твердой фазе. Однако ввиду отсутствия экспериментальных данных по теплотам сублимации твердого НДМГ для этого диапазона температур мы ограничились экстраполяцией параметров жидкой фазы на эти температуры.
Если учитывать атмосферную влажность, то каждая падающая капля будет представлять собой водный раствор НДМГ. Изменение массы капли и процентного содержания воды в ней происходит при падении капли за счет процессов испарения и конденсации молекул НДМГ и воды на поверхности капли. При этом формула (5) записывается отдельно для каждого компонента смеси.
Как и в предыдущих работах [1-3], предполагается, что температура поверхности капли Тв является квазиравновесной. Она определяется внутренними итерациями на каждом временном шаге при численном решении системы динамических уравнений (1)-(5) из условия баланса тепловой энергии, идущей на испарение молекул, и притока тепла к капле при ее аэродинамическом торможении и теплообмене с атмосферой за счет теплопроводности. Тепловая энергия, получаемая каплей при аэродинамическом торможении, принимается равной половине работы, совершаемой силой сопротивления воздуха. Этот механизм подвода тепла является определяющим на больших высотах при г > 80 км, в то время как теплопроводность наиболее существенна для малых высот, когда капли переходят в режим падения с относительно небольшими установившимися квазиравновесными скоростями. Обсуждаемые ниже результаты получены с использованием линейного интерполирования между двумя этими механизмами, аналогичного интерполированию в формуле (5).
Сила сопротивления воздуха, входящая в уравнения (I) ( I). также имеет две предельные формы. Для больших высот она квадратично зависит от скорости капли V относительно воздуха:
и = ~^кжг2крулг,
где р — плотность воздуха. В проведенных расчетах к = 1. Для малых скоростей на относительно небольших высотах применяется выражение для силы сопротивления, записанное в виде закона Стокса [1, 4]:
и = -бтг
Г в
Здесь ^ — динамическая вязкость воздуха; ^/^в — эмпирический поправочный множитель [4]. В расчетах для сшивки этих предельных случаев мы также использовали линейную интерполяцию.
Описанная физическая модель для обмена импульсом, энергией и массой между каплей и окружающей атмосферой на больших высотах и в переходной области высот, реализованная в текущем варианте программы, довольно упрощена, является предварительной и может быть уточнена путем привлечения имеющихся в литературе данных.
Для высот меньше 30 км и параметров атмосферы, таких как давление, температура, компоненты скорости ветра, коэффициенты турбулентной диффузии, влажность, использованы данные Гидрометцентра России, а для высот больше 30 км для давления и температуры использованы параметры Международной стандартной атмосферы (ГОСТ 4401-81), скорость ветра, влажность и коэффициенты диффузии полагаются нулевыми.
Разработанная программа позволяет рассматривать разовые выбросы топлива в любой точке траектории падения ступени начиная с момента отделения второй ступени. В основном режиме работы программы рассматривается последовательный выброс равных количеств топлива через равные промежутки времени. Тем самым в дискретной аппроксимации моделируется равномерное вытекание топлива из падающей ступени.
Параметры траектории падения капли НДМГ с начальным радиусом гк = 3 мм, образовавшейся при отделении второй ступени ракеты-носителя "Протон-М"
№ с г, км ж, км Г к, ММ т/то, % V, м/с г>г, м/с
1 0 127.790 0 3 100 4114.3 635.0
2 105.04 161.098 446.44 2.985 99.68 4341.5 0
3 244.22 102.603 1039.33 2.975 95.30 4452.2 -840.6
4 245.00 101.945 1042.69 2.971 94.84 4452.1 -845.1
5 252.56 95.398 1075.16 2.885 86.02 4436.7 -885.8
6 261.19 87.611 1111.74 2.500 55.25 4319.6 -912.2
7 269.75 80.000 1145.48 1.310 7.89 3701.3 -830.1
8 277.08 75.458 1164.29 0.547 0.59 1628.5 -402.4
9 293.19 71.919 1175.22 0.457 0.35 285.9 -137.0
10 329.95 68.000 1177.69 0.435 0.30 89.0 -88.8
И 504.43 61.213 1177.74 0 0 0 0
№ Кп М Не Т81° С Т, °С р, кг/м3
1 837.5 11.01 0.019 -104.9 127.2 1.62-10-«
2 10336 8.41 0.001 -115.5 390.6 1.32-10"9
3 21.21 16.09 1.339 -84.7 -82.7 6.44-10"7
4 18.93 16.15 1.509 -84.0 -84.0 7.22-10"7
5 5.927 16.64 5.091 -75.8 -96.1 2.38-10"6
6 1.712 16.20 17.16 -65.4 -96.1 9.49-10"6
7 0.926 13.45 25.75 -58.3 -84.8 3.35-10"5
8 1.099 5.791 9.193 -76.4 -76.4 6.76-10"5
9 0.775 1 2.223 -99.1 -69.8 1.15-10"4
10 0.471 0.304 1.093 -97.3 -59.7 1.98-10"4
И оо 0 0 -84.8 -42.1 4.83-10"4
Приведем результаты расчета падения капли НДМГ с начальным радиусом 3 мм, образовавшейся сразу после отделения ступени на высоте г = 127.790 км. Результаты, соответствующие нескольким ключевым моментам движения капли, представлены в отдельных строчках таблицы. В таблице приведены следующие параметры: £ — время в секундах с момента отделения ступени; г — высота капли над уровнем моря; х — удаление капли от места отделения ступени по дуге на уровне моря; гк — радиус капли; т/т0 — масса капли в процентах к начальной массе; V — скорость капли; уг — вертикальная компонента скорости капли; Кп — число Кнудсена; М — число Маха; Ее — число Рейнольдса; Т — температура поверхности капли; Т — локальная температура воздуха; р — локальная плотность воздуха.
После выброса жидкого топлива его поверхность охлаждается за счет испарения до квазиравновесной температуры Т3. Предполагаем, что начальная капля радиуса гк = 3 мм уже имеет установившуюся квазиравновесную температуру Т = -104.9 °С (1-я строка в таблице) и соответствующую массу 0.104 г. Начальная скорость капли V = 4414.3 м/с совпадает со скоростью ступени. Далее капля совершает баллистический полет в условиях крайне разреженной атмосферы, поднимаясь через время £ = 105.04 с до максимальной высоты своей траектории г = 161.098 км (2-я строка). При этом капля охлаждается еще па 10.6°. Потеря массы капли па этот момент времени за счет испарения составляет всего лишь 0.32 % от начальной массы, и капля движется горизонтально со скоростью 4341.5 м/с. Удаление этой точки от начальной точки по дуге на уровне моря в системе отсчета Земли составляет х = 446.44 км.
При дальнейшем полете капли в баллистическом режиме ее высота начинает уменьшаться и полная скорость возрастает. Но возрастает также и сила торможения со стороны воздуха, поэтому в момент времени £ = 244.22 с на высоте г = 102.603 км скорость капли достигает своего максимального значения — 4452.2 м/с. Это соответствует числу Маха М = 16.09. Вертикальная компонента скорости — уг = -840.6 м/с. Потеря массы НДМГ составляет уже 4.70 %, а температура поверхности капли Т3 = -84.7 °С, что всего лишь на 2.0° ниже температуры окружающего воздуха па этой высоте.
Очевидно, что полученные численные результаты практически не зависят от того, в какой момент времени произведен выброс капли из ракетной ступени на рассмотренной части баллистической траектории, поскольку речь идет о движении по инерции в очень разреженной атмосфере. То же самое относится и к каплям других начальных размеров.
В процессе дальнейшего движения температура капли продолжает повышаться. На высоте г = 101.945 км (4-я строка) температура поверхности капли сравнивается с локальной температурой атмосферы Т3 = Т = -84.0 °С. Эту точку верхней части термосферы условно можно считать началом активной фазы торможения капли. Число Маха растет и достигает своего максимального значения М = 16.64 при высоте г = 95.398 капли.
При продолжающемся уменьшении полной скорости вертикальная компонента скорости все еще увеличивается под действием силы тяжести и достигает своего максимального значения 912.2 м/с на высоте г = 87.611 км (6-я строка). На этот момент времени осталось лишь 55.25 % от начальной массы капли.
Аэродинамический нагрев капли при усиливающемся торможении ведет к росту тем-
-58.3 °
26.5°
связано также с отсутствием активного теплообмена между каплей и атмосферой за счет теплопроводности. Нагрев капли сопровождается усилением ее испарения, поэтому до высоты 80 км долетает только 7.89 % от начальной массы капли и, соответственно, радиус капли уменьшается более чем в два раза — до значения гк = 1.310 мм.
Уменьшение скорости капли ведет к уменьшению скорости подвода тепловой энергии к капле, что вызывает последующее понижение ее температуры. На высоте г = 75.458 км (8-я строка) температура поверхности капли снова сравнивается с локальной температурой атмосферы Т8 = Т = -76.4 °С. Здесь полная скорость капли уже уменьшилась почти в три раза — до величины 1628.5 м/с при числе Маха М = 5.79. Эту точку мезосферы условно можно считать завершением фазы наиболее активного аэродинамического торможения капли. Масса капли составляет 0.65 % от исходной массы, и радиус капли уменьшился почти в шесть раз — до гк = 0.547 мм.
Еще через 16.11 с полета на высоте г = 71.919 км (9-я строка) происходит переход через звуковой барьер при скорости 285.9 м/с. При этом радиус капли уменьшается до значения 0.457 мм. Температура капли Т8 = -99.1 °С, что уже па 29.3° ниже температуры окружающей атмосферы.
На высоте 68.000 км (£ = 329.95 с, 10-я строка) в средней части мезосферы горизонтальная составляющая скорости уже почти исчезает и капля падает практически вертикально со скоростью 89.0 м/с. Устанавливается режим падения с квазиравновесной
скоростью. Радиус капли при этом равен 0.435 мм. Температура капли Т8 = -97.3 °С,
°
Окончательное испарение капли происходит примерно через три минуты на высоте 61.213 км (£ = 504.43 с, 11-я строка) при температуре капли Т8 = -84.8 °С. Удаление этой точки от места образования капли по дуге на уровне моря вдоль земной поверхности составляет 1177.74 км.
Представленные результаты иллюстрируют тот факт, что выброс остатков топлива из бака ракетной ступени на больших высотах, когда ступень имеет еще достаточно большую скорость, решает проблему устранения экологических рисков, связанных с возможностью загрязнения земной поверхности и приземного слоя атмосферы остатками ракетного топлива.
Помимо загрязнения земной поверхности топливом, долетевшим до земли в капельном виде, интерес представляет также моделирование временной эволюции облака паров топлива, образовавшегося при испарении в атмосфере падающих капель. Динамика этого облака определяется ветровым сносом и атмосферной турбулентной диффузией и описывается трехмерным уравнением диффузии с переносом.
Задача численного моделирования временной эволюции в атмосфере облака пассивной примеси сама по себе не нова, но специфику данной задаче придает наличие подвижного локализованного источника капель (падающая ступень ракеты) и подвижного распределенного источника паров топлива (падающие в атмосфере испаряющиеся капли). Это приводит к тому, что при использовании разностных методов фиксированная в пространстве сетка уже не дает возможности оптимально описывать рассматриваемые физические процессы. Поэтому можно ожидать, что использование подвижных адаптивных сеток может заметно улучшить точность моделирования таких процессов.
Нами рассмотрена возможность применения к данной задаче универсального метода построения адаптивных сеток, описанного в [5, 6] и основанного на численном решении обращенных уравнений Бельтрами и диффузии.
При применении такого метода построения сеток осуществляется переход к новым пространственным переменным (С2,£3), в которых уравнение переноса и диффузии аппроксимируется разностным уравнением на простой кубической пространственной сетке. В дальнейшем временная задача решается в переменных 2,£3). Узлам ку-
бической сетки в переменных (С^^3) соответствуют узлы криволинейной сетки в физических переменных (х^х^х3). Вектор-функция х(£) находится итерационно на каждом временном шаге как решение обращенного сеточного уравнения диффузии [5, 6]. Использование в сеточном уравнении подходящей мониторной функции, зависящей от плотности примеси, позволяет получить необходимые сгущения сетки в физическом пространстве в областях с повышенной плотностью примеси.
Эволюция во времени плотности пассивной примеси р(£, х) определяется дифференциальным уравнением, которое в физических пространственных переменных (ж1, х2, ж3) в дивергентной форме имеет вид
где иг(Ь, х) — компоненты вектора скороети ветра; В^ (£, х) — матрица коэффициентов турбулентной атмосферной диффузии; Г(£, х) — плотность источника примеси. Здесь и далее, как обычно, по повторяющимся индексам предполагается суммирование. Для уравнения (6) задаются начальная плотность р(0, х) и некоторые краевые условия на границе рассматриваемой пространственной области.
(6)
При переходе к новым пространственным переменным (£ 1,С2,С3) уравнение (6) принимает форму
д(1р) + д_ Г д?
д£ дСЧ дхк
где I — якобиан преобразования.
дь ) кздедх>
дхг
/ = ёе! ( ——
1^, (7)
де
Авторами работ [5, 6] для нахождения эффективного отображения х ^ е развит подход, основанный на методах римановой геометрии. В данной задаче мы рассматриваем этот подход в следующем варианте.
Предполагается, что значения вектор-функции х(£) на границе рассматриваемой области заданы и в каждый момент времени £ функция х(£) является решением краевой задачи для обращенного уравнения диффузии, имеющего вид
ЗУ дю(р)дСп= ( )
ю[р) дхк дхк д£рд£1 д£т дхг { 1
Здесь -ш(р) — эмпирическая сферическая мониторная метрика, подбором которой можно добиться желательного сгущения (или разрежения) сеточных узлов в физическом пространстве в областях больших (или малых) значений плотности р(£, х).
Задача решается дискретными шагами по времени. Рассмотрим один такой временной шаг величины Д£. Пусть началу временного шага соответствует значение времени ¿0) тогда концу временного шага отвечает момент времени ¿0 + Д£-
Решение полной задачи (6)-(8) па временном шаге Д£ разбивается па два этапа. Вначале проводится адаптация сетки, т. е. ищется функция х(£) как решение краевой задачи для уравнения (8). При этом плотность в физическом пространстве считается фиксированной р(х) = р(£0, х). На втором этапе на полученной сетке решается уравнение (7) и находится плотность р(£0 + Д£, х). На этом временной шаг заканчивается, и переходим к следующему временному шагу.
На первом этапе временного шага для нахождения решения эллиптического уравнения (8) используется метод установления, т. е. вместо уравнения (8) рассматривается его параболизованная форма:
дх% _ 2
=
, л дер де? д2хг дш(р) дег
и)(р)-
дхк дхк деРде« дет дх^
(9)
Здесь через т обозначен итерационный параметр, по которому и происходит установление.
е
аппроксимируются центральными разностями. При этом для расчета элементов обратной матрицы Якоби, входящих в (9), используется выражение
д£р _ 1 ~ 1
дхк+1 дхк+2 дхк+1 дхк+2
дер+1 дер+2 дер+2 дер+
(10)
где в правой части значения верхних индексов понимаются в циклическом смысле, т. е. 4 заменяется на 1, а 5 — на 2.
р
Поскольку построение новой сетки на очередном временном шаге является лишь вспомогательной задачей для решения основной задачи описания временной эволюции распространения примеси, то точность нахождения "стационарного" решения для сеточного уравнения (9) может быть относительно небольшой. Соответственно, число итераций будет также небольшое.
Но одной только функции х(£), полученной па кубической сетке при решении уравнения (9), еще недостаточно, чтобы перейти к нахождению решения уравнения (7). Требуется определить изменившиеся значения плотности р(£0, £) в узлах кубической сетки. Для этого одновременно с уравнением (9) на первом этапе временного шага находится эта функция р(£0, £) как результат решения уравнения
д(1р) д
дт д?
1р
д£г дх
0, (П)
дхк дт
которое эквивалентно условию независимости плотности р(т, x) от итерационного пат
= 0. (12)
дт
Итак, задача полностью решается на простой кубической сетке в ^-пространстве. При этом на первом этапе временного шага проводятся внутренние итерации (по парат
для уравнения (7).
Изложенный подход к использованию адаптивных сеток сформулирован в результате подробного обсуждения с В. Д. Лисейкиным, И. А. Васевой и Ю. В. Лихановой. Автор выражает также признательность А. И. Бородулину и Б. М. Десяткову (НИИ аэробиологии ГНЦ ВБ "Вектор") за расчет параметров траектории движения отделившейся второй ступени ракеты-носителя "Протон-М" и полезные обсуждения рассматриваемых в данной работе проблем.
Список литературы
flj Климова Е.Г., Мороков Ю.Н., Ривин P.C. и др. Моделирование загрязнения поверхности земли ракетным топливом // Оптика атмосферы и океана. 2004. Т. 17, № 9. С. 769-773.
[2J Климова Е.Р., Мороков Ю.Н., Ривин P.C. и др. Математическая оценка зон загрязнения поверхности земли ракетным топливом при падении отделяющихся частей ракет-носителей // Оптика атмосферы и океана. 2005. Т. 18, № 5-6. С. 525-529.
[3] Шокин Ю.И., Ривин P.C., Климова Е.Г., Мороков Ю.Н. Математическое моделирование загрязнения территории Восточного Казахстана и Алтая ракетным топливом при падении отделяющихся частей ракет-носителей // Вестник КазНУ. Сер. Математика, механика, информатика. 2005. Т. 18. С. 93-102.
[4J Александров Э.Л. О поведении капель ракетного топлива в атмосфере // Метеорология и гидрология. 1993. № 4. С. 36-45.
[5J Методы римановой геометрии в задачах построения разностных сеток / Ю.И. Шокин,
В.Д. Лисейкин, A.C. Лебедев и др. Новосибирск: Наука, 2005. 256 с. [6J Построение разностных сеток с помощью уравнений Бельтрами и диффузии / А.Г. Глас-сер, В.Д. Лисейкин, Ю.И. Шокин и др. Новосибирск: Наука, 2006. 184 с.
Поступила в редакцию 14 марта 2008 г.
