CC BY
59
Вычислительные технологии
Том 8, № 2, 2003
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИЧЕСКИХ ПРОЦЕССОВ В СОЛНЕЧНОЙ АТМОСФЕРЕ
Д. В. РОМАНОВ, К. В. РОМАНОВ Красноярский государственный торгово-экономический институт,
Россия e-mail: ds@kgtei.kts.ru
The article is devoted to numerical modelling of non-stationary processes in the solar atmosphere. For the description of atmospheric plasma the choice of approximation of one-fluid gas dynamics taking into account the radiation heat conductivity and thermal conduction processes is proved. The explicit conservative scheme for a numerical solution of the written out equations set is constructed. The used interpolation algorithm of volumetric forces allows to reduce errors of account stationary solutions found by relaxation method.
The numerical results of test problem and the solar atmospheric heating with overturning acoustic waves problem are indicated.
Введение
В настоящее время при изучении структуры солнечной атмосферы необходимо использовать методы численного моделирования, так как невозможно аналитически учесть все значимые нелинейные физические процессы. В качестве примера можно привести следующие задачи: расчет распределения по высоте газодинамических параметров атмосферы решением обратной задачи по данным наблюдений в различных спектральных диапазонах [1]; прогрев солнечной атмосферы опрокидывающимися при распространении в стратифицированной среде акустическими волнами [2-4]; отклик атмосферы на импульсное термическое или динамическое возмущение [5, 6]; эволюция сложных магнитных структур в активных областях солнечной атмосферы [7] и многие другие [8].
Настоящая работа посвящена особенностям моделирования газодинамических течений в солнечной атмосфере. Ранее подобные задачи решались в одномерной постановке, недостаточной для адекватного описания нестационарных физических процессов, в значительной степени неоднородных по всем направлениям [8]. В настоящей работе построена разностная схема для решения двумерной системы уравнений газовой динамики на базе двумерной схемы [9] (используемой для решения уравнений мелкой воды) с добавлением искусственной вязкости с учетом процессов кинетического и лучистого теплопереноса и измененной интерполяции гравитационных сил.
В разд. 1 приведены значения физических параметров солнечной атмосферы и обоснован выбор математической модели для ее описания. В разд. 2 выписана двумерная
© Д. В. Романов, К. В. Романов, 2003.
система уравнений газовой динамики с учетом вязкости, теплопроводности, гравитации и объемных лучистых потерь. Решению полученной системы уравнений методом конечных разностей посвящены разд. 3-5. В разд. 6 исследуется вопрос о невязке стационарных решений, полученных методом установления. Результаты тестовых расчетов приведены в разд. 7. В качестве примера использования полученной схемы в работе исследуется задача о прогреве солнечной хромосферы и короны опрокидывающимися акустическими волнами. Результаты расчетов представлены в разд. 8. В заключении дан обзор полученных результатов.
1. Физические параметры солнечной атмосферы
Солнечная атмосфера состоит в основном из водородно-гелиевой плазмы, химический состав которой приведен в табл. 1. Относительная концентрация элементов в солнечной фотосфере нормирована так, что п(Я1) = 106 см-3 [10]. Элементы с концентрацией меньше 104 в таблицу не включены. Диапазон изменения физических параметров солнечной атмосферы с ростом высоты велик. Это следует из табл. 2, где приведены характерные значения параметров плазмы атмосферы в зависимости от высоты.
Атмосферу Солнца принято делить по свойствам на три области. Основанием атмосферы является фотосфера. Ее нижней границей служит видимая поверхность Солнца с температурой 6600 К, а верхней границей — зона температурного минимума, где температура падает до 4140 К по данным модели [1]. Выше идет хромосфера, в которой температура возрастает до значений порядка 105 К. Область над хромосферой названа короной. Все эти области сильно различаются по своим физическим характеристикам. Общим качеством является очень маленькая оптическая плотность среды, делающая пригодным описание лучистого переноса в рамках оптически тонкой среды (объемные потери энергии). Фотосфера — наиболее плотная и холодная часть атмосферы, она состоит практически из неионизованной плазмы. Теплопроводность и вязкость в ней пренебрежимо малы (см. табл. 2). Плазма хромосферы уже заметно горячее, и коэффициент ионизации достаточен для того, чтобы увеличить коэффициент теплопроводности до величины, достаточной для влияния на профиль температуры в верхней части хромосферы, где плотность плазмы очень мала. Излучение из хромосферы идет в основном в ультрафиолетовом диапазоне. В короне теплопроводность играет доминирующую роль. Данная модель предназначена для исследования подкорональной части атмосферы Солнца из-за явного описания теплопроводности.
Основные параметры солнечной атмосферы приведены в табл. 2. Здесь ^ — молярная
Таблица1
Химический состав водородно-гелиевой плазмы
Элемент 3 п, см 3 Элемент -3 п, см 3
И 3.32 ■ 1010 Mg 9° 5 1 о Сл
Ив 1.96 ■ 109 А1 3 1 О
С 1.17 ■ 107 6 0 1 .0 1
N 3.16 ■ 106 Я 6.17 ■ 105
0 2.04 ■ 107 Са 7.08 ■ 104
N6 2.46 ■ 106 Рв 2.19 ■ 105
сл 0 1 о N1 1.95 ■ 104
Т аблица2
Характерные параметры плазмы солнечной атмосферы на различных высотах
Л-, см р, эрг/см3 т, к p, г /см3 ^, г/моль 7 д, см /с2
2.00 • 1010 1.10 • 10-3 2.60 106 3.13 10-18 0.617 1.63 1.65 • 104
4.00 • 109 3.50 • 10-3 1.85 106 1.41 10-17 0.620 1.63 2.45 • 104
1.00 • 109 5.02 • 10-3 1.31 106 2.93 10-17 0.635 1.63 2.66 • 104
4.05 • 108 5.64 • 10-3 9.30 105 4.75 10-17 0.651 1.63 2.71 • 104
2.64 • 108 5.90 • 10-3 5.11 105 9.42 10-17 0.678 1.63 2.72 • 104
2.48 108 6.14 • 10-3 5.91 104 1.05 10-15 0.838 1.63 2.72 • 104
2.29 108 7.96 • 10-3 2.00 104 4.30 10-15 0.897 1.63 2.72 • 104
1.85 108 6.52 • 10-2 7.80 103 1.19 10-13 1.18 1.63 2.73 • 104
8.63 • 107 6.32 101 5.39 103 1.82 10-10 1.29 1.63 2.73 • 104
8.01 • 107 1.05 102 5.15 103 3.19 10-10 1.29 1.63 2.73 • 104
5.44 • 107 1.21 103 4.14 103 4.55 • 10-9 1.30 1.63 2.74 • 104
1.97 • 107 3.26 104 5.06 103 1.00 • 10-7 1.30 1.63 2.74 • 104
3.03 106 1.14 105 6.66 103 2.60 • 10-7 1.26 1.58 2.74 • 104
Л-, см Л, см а А, см Ие С1
2.00 • 1010 2.12 1010 1.00 6.78 • 108 1.24 3.17 • 10-2
4.00 • 109 1.01 1010 1.00 8.72 • 107 4.60 1.34 • 10-1
1.00 • 109 6.43 109 1.00 2.27 107 1.11 • 101 3.42 • 10-1
4.05 108 4.38 109 1.00 7.37 106 2.30 • 101 7.20 • 10-1
2.64 108 2.30 109 1.00 1.16 106 7.55 • 101 2.34 10
2.48 108 2.16 • 108 1.00 1.41 103 5.21 103 1.32 102
2.29 108 6.80 • 107 8.66 • 10- 1 4.21 101 5.32 • 104 5.17 102
1.85 108 2.01 • 107 2.23 • 10' 5 2.77 10 2.08 105 5.68 102
8.63 • 107 1.27 107 1.70 • 10“ 9 1.81 • 10-3 1.92 • 108 4.81 105
8.01 107 1.21 107 4.06 • 10- 10 1.03 • 10-3 3.21 • 108 8.01 105
5.44 107 9.72 106 1.87 • 10- 13 7.24 • 10-5 3.68 109 9.20 106
1.97 107 1.19 107 2.32 • 10- 10 3.28 • 10-6 9.90 1010 2.47 108
3.03 106 1.60 107 5.88 • 10“ 7 1.27 • 10-6 3.48 1011 8.71 108
масса; 7 — показатель адиабаты; Л = КТ/ цд — шкала высот [8] (высоты отсчитываются от радиуса г0 = 6.959 • 1010 см); а = Пг/иа — степень ионизации плазмы; Л — характерная длина пробега частицы; И,е — число Рейнольдса; С1 — число Клаузиуса. Таблица построена по данным работ [10-12]. Этих данных достаточно, чтобы выбрать математическую модель для описания физических процессов, протекающих в солнечной атмосфере. Во всем диапазоне высот, представленном в табл. 2, длина пробега частиц Л много меньше характерного пространственного масштаба — шкалы высот Л1, поэтому применимо приближение газовой динамики. Вязкостью в фотосфере и хромосфере (Т < 105 К) можно полностью пренебречь, а теплопроводность следует учитывать только в областях с достаточно разреженной плазмой: коэффициент теплопроводности не зависит от плотности, в то время как энергосодержание в единице объема ей пропорционально, из-за чего роль теплопроводности при фиксированной температуре растет с уменьшением плотности плазмы.
В атмосфере Солнца большое значение имеют процессы лучистого теплопереноса. Длина пробега излучения, как правило, очень велика, и характерное значение оптической
ХЛ = р/(рд) = р/(йр/йг) — диапазон высот, на котором плотность и давление уменьшаются в е раз из-за действия силы тяжести в условиях гидростатики.
толщины в ней много меньше единицы (оптическая толщина набирает значение около единицы в тонком слое толщиной порядка 100 км вблизи фотосферы [4]). Следовательно, для описания процессов лучистого переноса в первом приближении можно использовать приближение объемных потерь [11-13]. При описании лучистого теплопереноса без точного решения уравнений переноса для всего спектра (что повышает размерность задачи на два) в переходной зоне (к £ [0.500] км) можно использовать эддингтоновское приближение, корректно описывающее предельные случаи оптически тонкой и толстой сред [14].
Солнечная атмосфера не находится в состоянии локального термодинамического равновесия (ЛТР), поэтому степень ионизации плазмы, коэффициенты поглощения и излучения света необходимо определять непосредственно из уравнений баланса ионизации и рекомбинации с учетом всех значимых микроскопических процессов без использования принципа детального равновесия. В атмосфере Солнца ионизация осуществляется в основном электронным ударом, а рекомбинация реализуется за счет фоторекомбинации вместо обратного процесса тройной рекомбинации. Это равновесие носит название коронального, и в этом случае степень ионизации плазмы определяется по формуле Эльверта [15]. Мощность объемных потерь для случая коронального равновесия найдена в работах [11-13]. Средние по спектру (росселандовы) коэффициенты поглощения в литературе приведены только для случая ЛТР [16, см. библиографию], поэтому их использование затруднено тем обстоятельством, что коэффициенты поглощения и излучения уже не связаны соотношениями Кирхгофа [17].
По наблюдательным данным, солнечная атмосфера является в значительной степени нестационарным объектом. Основные источники возмущений — конвективные течения (наблюдаемые как грануляция), вспышечные процессы и магнитные поля (солнечные пятна, вспышки, корональные петли), а также выбросы вещества в атмосферу (протуберанцы, корональные транзиенты, спикулы) [8]. Все эти явления в значительной степени неоднородны, и для их адекватного описания желательно решать задачу как минимум в двумерной постановке.
2. Описание системы уравнений
Выпишем систему уравнений одножидкостной газодинамики с учетом вязкости, теплопроводности, действия сил гравитации и объемных лучистых потерь:
Здесь U — вектор газодинамических параметров среды; F — вектор потока вдоль оси х; G — вектор потока вдоль оси y; W — вектор, описывающий объемные потери; L — вектор, описывающий действие гравитации. Индекс id отвечает недиссипативной компоненте потоков, visc — вязким компонентам потоков, heat - слагаемым, соответствующим теплопроводности. Покомпонентная запись этих векторов следующая:
(1)
глт С 2 { р(и2 + V2)
= 4 р^; риь; р + р^ ; V I р + ре +---------^----
т Г„ ^ ди ъ {2 Л п. {ди дv
г- = \0; -2"{Х Чз" - С; Лт У; -ЧдУ + дх
ди2 {2 " {ди дv
-пГХ + Ч зп - ^Л1Т " - + дХ
,т Г„ {ди , дV" ^ д" , {2
еу* = |°; -п {ду + дХ] ; -2пдУ + (,зп - С 1 Шу У;
дV2 {2 " {ди д"
-г% + Ч 3п - V Лу У - ^ 1ду + дХ
^т Г д^ _т Г дТ
Fheat = |0; 0; 0; -Xа^J ’ ^Ьеа1 = |0; 0; 0; -Хду
Ьт = {0; рдх; рду; р(идх + vду)} ,
Wт = {0; 0; 0; -Жгаа(Т,р)} .
В этих уравнениях р, р, е, Т обозначают плотность, давление, удельную (на единицу массы) внутреннюю энергию и температуру плазмы; V = (и, V) — скорость среды; п, С — коэффициенты вязкости [18]; х — коэффициент теплопроводности; g = (дх, ду) — вектор ускорения свободного падения; Жгаа — мощность объемных лучистых потерь. В случае коронального равновесия и х, и Жгаа являются функциями только от температуры.
3. Разностная схема
Большой диапазон изменения параметров среды (ртах/ртт ^ 1022 для равновесной изотермической атмосферы высотой порядка 104 км при Т ~ 4000 К) накладывает ограничения на алгоритм, используемый при обращении матрицы разностной системы уравнений, особенно с использованием неявных методов. Для сохранения точности алгоритм приходится переписывать таким образом, чтобы исключить вычисления отношения двух отличающихся на малые значения величин, суммирования больших и малых величин и т.д. (см., например, [19]). С учетом этого факта и того, что вызванные наличием гравитации и лучистого теплопереноса дополнительные ограничения на шаг по времени бывают более сильными, чем критерий Куранта, для решения динамической части системы уравнений лучше всего выбрать явную консервативную схему.
Из-за наличия гравитации среда является стратифицированной, поэтому даже слабые возмущения при распространении вверх имеют амплитуду, растущую по экспоненциальному закону. Как следствие этого, используемая схема должна хорошо описывать ударные волны. Для этого существует несколько подходов: использование схем с большой счетной вязкостью (схемы Лакса — Фридрихса), использование ТУБ-модификаций схем с большим порядком аппроксимации и использование искусственной вязкости для сглаживания осцилляций на фронтах ударных волн.
Использование композитных схем [9] и схем типа Лакса — Фридрихса затруднено тем обстоятельством, что в схеме Лакса — Фридрихса происходит двойное усреднение:
ттга+1/2 = 1 (\тп + ип + ип + ип ^ +
иг+1/2,^+1/2 = 4 \иг^ + иг+1,^ + иг,^+1 + иг+1,^+^ +
-П+1/2 . и уп+1/2 . ууП+1/2 . ууП+1/2
І-1/2.-1/2 + иг+1/2,і-1/2 + иг-1/2.+1/2 + иг+1/2.+1/2
Из-за этого на стационарном решении, где р, р ж ехр(—г/Л), происходит нефизический рост массы со скоростью, которую нетрудно оценить:
Этот эффект обусловлен исключительно кривизной стационарного решения, и его наличие исключает использование схемы типа Лакса — Фридрихса для решения системы уравнений вида (1) без предварительной доработки. Одним из методов коррекции является интерполирование в полуцелые точки с учетом гравитации
но это не устранит недостаток полностью, только уменьшит скорость роста массы до величины, пропорциональной ііА/ііх. Так как для применения композитных схем необходимо сначала создать быструю монотонную консервативную схему с хорошей счетной вязкостью, но без счетного роста массы, использование искусственной вязкости для построения неосциллирующей на ударных волнах схемы является самым простым и надежным методом. В связи с этим для решения динамической части системы уравнений выбрана явная схема второго порядка аппроксимации, переходящая в одномерном случае в схему Лакса
— Вендроффа. Подобная схема используется для подшага в композитной схеме, предназначенной для решения уравнений мелкой воды в работе [9].
Основная идея построения схемы базируется на результате работы [20], где показано, что при использовании метода Годунова в двумерных задачах для определения средних по ячейке значений величин на сдвинутой сетке достаточно решить одномерные задачи на границах ячейки. Опишем шаг по времени с учетом одной гравитации (подробнее об аппроксимации диссипативных потоков и лучистого переноса написано в следующих разделах).
Шаг по времени разбит на три этапа. На первом находят значения и
из решения одномерной задачи на границе ячейки:
п+1/2 _ 1 + ехр(±^ж/Л)
РІ±1/2 _ РІ 2
п+1/2 _ РІ±1/2 _ Рі
Рі(1 + 0.5еЬ (кх/Л)),
На втором шаге определяются значения ип+11//22 .
і+1/2, .+1/2:
+ 2 Ч1 (Ц,* + и?+1-* + и«+‘ + и?+^+1
И заключительный шаг делается следующим образом:
И*1 = И ■+
г,* * 1
+ 5Г ^<И+1А;+1/2> + F(И—Х—1/2) - ^(И+%+1/2) - ^«й*—1/2)) +
2Ъ.У (С(ИГ-+Й,-1/2) + С(И”++11//2*-1/2) - С(И“-Х+1/2) - С(И”++11//2*+1/2^ +
+Л( 1 (Иг,* + И*ч) •
(2)
Формально это уравнение неявное из-за члена Е((И™*- + И+1)/2), но его значение можно найти прямыми вычислениями, так как первая компонента вектора Ь равна нулю, а остальные зависят только от ранее найденных компонент вектора И.
4. Учет вязкости
Как следует из данных табл. 2, значение вязкости в солнечной атмосфере очень мало, за исключением области верхней короны. Тем не менее для сглаживания осцилляций разностного решения около разрывов в схему добавляется небольшая искусственная вязкость. Вязкие слагаемые системы уравнений можно разделить на две части — одна часть содержит производные параметров плазмы только по направлению потока, что дает после вычисления дивергенции производные первого и второго порядка только по одной переменной (индекс 1), в то время как вторая группа слагаемых содержит оставшиеся члены, дающие при определении дивергенции потока смешанные производные (индекс 2):
гТ:
сТ
г:
Т
.4 \ ди ду ду /4 \ ди
0; - -п + С ; -пт-; -пу^- - + С
3 дх дх дх \ 3 дх
ди /4 Л ду ди /4 Л ду
0; -,'ду;4 3п + Чту; -п% - Ч 3п + С]ду
2 ду ди 2 ду ди
0; ип - С)зу; -пду; Ч3п - С)дУ - пуду
ду 2
'дХ; 13
0; -пдх; [-п - ; у I 3п - С) — ~'пи
ди
дх
ду
дх
Входящие в потоки Г2, С2 производные аппроксимируются через И™+11//22*±1/2:
ди п+1/2 п+1/2 иг+1/2, *+1/2 - иг—1/2, *+1/2 ду п+1/2 п+1/2 7 7 ' 71 ' уг+1/2,*+1/2 уг—1/2, *+1/2
дх ди Ъ ’ г, *+1/2 Ъх п+1/2 п+1/2 7/ 7 ?/ иг+1/2,*+1/2 иг+1/2, * — 1/2 дх ду Ъ г, *+1/2 уга+1/2 уп+1/2 уг+1/2,*+1/2 уг+1/2,* —1/2
ду Ъ ’ г+1/2, * ЪУ дх Ъ г+1/2, * ЪУ
2
Входящие в потоки Г і, Сі производные аппроксимируются через и™?., и™±1 3, ип3-±1-
ди ип _ ип иі+1,3 иі,3 дь
дх і+1/2,3 кх дх
ди ип _ ип і, 3+1 - і, 3 дь
дУ і, 3+1/2 НУ ’ дх
і+1/2,3
1)П ________
ьі+1,3 і, 3
кх
Ьп ___ 7?п
ьі,3+1 ь
і, З
ку
г, 7+1/2 '°У
Все производные (в том числе и разностная дивиргенция потоков Г1, ..., 02) вычисляются центральными разностями, но потоки со смешанными производными находятся на полуцелом, а потоки со вторыми — на целом слое по времени. Консервативность системы от этого не страдает, а уменьшение порядка аппроксимации по времени до первого несущественно, так как искусственная вязкость вводится только для уменьшения осцилляций решения за фронтами ударных волн.
В процессе счета вычисление диссипативных потоков и их дивергенции происходит на
г±1
заключительном шаге (3), после определения и"^1/2 з-±1/2.
5. Лучистый теплообмен и теплопроводность
В настоящей работе в основном исследуются процессы в оптически тонких частях атмосферы Солнца — хромосфере и короне, поэтому процесс лучистого теплообмена хорошо описывается в приближении объемных потерь. Мощность излучения из единицы объема плазмы заданного в табл. 1 химического состава, найдена в работах [11-13].
В процессе численного интегрирования слагаемые, отвечающие теплопереносу, из-за теплопроводности и процессов излучения находятся явным образом. Обозначим компоненты вектора потока тепла как q = (я(х), 3(У)), тогда
Пт/2,3 = -х (0-5(Тт,з + Ті,з)) 1+1,3 і3
^+1/2 = -х (0.5(Ті,з+1 + Ті,з)) ь
Ті,3+1 - Ті,3
ип+1 - ип ■
і, 3 і, З
т
П(х) - П(х) П(у) - П(у)
Уі+1/2,3 пі-1/2,3 Пі, 3+1/2 Пі,3-1/2
кх
ку
- ^гаа(Т)пнП.
Здесь w обозначает последнюю (четвертую) компоненту вектора и; пн = Рн/тн — концентрация водорода; пе — концентрация электронов. Концентрация электронов находится из условия баланса между ионизацией электронным ударом и фоторекомбинацией [12], ее зависимость от температуры приведена на рис. 1. Зависимость коэффициента теплопроводности от температуры дана на рис. 2. При небольших температурах (пе ^ пн) теплопроводность определяется нейтральными атомами, при больших температурах (пе ~ пн) основной вклад дают электроны. Зависимость мощности объемных потерь (приведенной к пепн) от температуры изображена на рис. 3.
6. Погрешности стационарного решения
Так как в схеме Лакса — Вендроффа значения на полуцелом слое по времени используются только для определения потоков, эта схема избавлена от описанного в разд. 3 недостатка —
Рис. 1. Зависимость концентрации электронов от температуры.
Рис. 2. Зависимость коэффициента теплопроводности от температуры.
Рис. 3. Зависимость мощности объемных потерь от температуры [12].
счетного роста массы. Тем не менее кривизна стационарного решения приводит к другому разностному эффекту — в разностном стационарном решении скорость отличается от нуля. Действительно, условием стационарности для динамической части системы уравнений является равенство нулю дивергенции потоков на полуцелом слое по времени:
(р (иГЙ,^) + р (И-Х-^) - пи^1/^) - р (И+Х-1/2)) +
2Л,х
+ ^ (С(И"_+!/27-1/2) + е(И”++11//2у_1/2) - С(И“+1/22, + 1/2) - С(И”++1/2у+1/2)) +
+т ь( 2 (и;:,+и“+ч) =о.
Поэтому скорость равна нулю на полуцелых слоях по времени, а на целых она будет отличаться от нуля на величину порядка О (т/2). Профиль скорости разностного стацио-
и и ТТ Т П+1/2
нарного решения при использовании линеинои интерполяции и при определении ^+1/2 приведен на рис. 4. Безразмерные параметры задачи следующие: д = -3, Нх = Ну = 10-2, т = 10-4, п = С = 2• 10-3 • р, р0 = 0.99985, р0 = 0.99991, а = 0.5. Видно, что скорость в целых узлах по времени равна малои постояннои величине, а ее отличие от нуля вызывает осцилляции решения на границе.
Этот эффект не имел бы большого значения, если бы найденные на целых слоях по времени значения скорости не использовались, внося соответствующую ошибку, для определения величин р, Т, е в уравнениях теплопереноса. С целью исправить этот недостаток была изменена процедура интерполяции при определении Ь в полуцелой точке.
Рис. 4. Разностное стационарное решение при фиксированной температуре Т(х) = 0.156347.
Рис. 5. Стационарное разностное решение при использовании интерполяции (4). Остальные параметры такие, как на рис. 4.
Сформулируем вспомогательную одномерную задачу: пусть при решении задачи на отыскание равновесного распределения параметров мы используем схему Лакса — Вен-дорффа с некоторой линейной интерполяцией при определении Ьг+1/2. Разностное стационарное решение должно удовлетворять следующим условиям: скорость равна нулю, температура постоянна, давление и плотность зависят от координаты по показательному закону [8]
Рг = РоК\ рг = РоХ\ (3)
Шаг по времени делается как обычно. С учетом того, что и™ = 0,
рп+1/2 = 1(рП . рП)
Рг+1/2 = 2 (Рг+1 + Рг )}
(ри)П+11/22 = 2кг (^(1 - к) + адро +(1 - а)д^ро)
Рис. 6. Стационарное разностное решение при фиксированном распределении температуры Т(х) = 0.156347 + 0.05(1 - шз(пх/2)), ро = 1.0068, ро = 1.0068, остальные параметры, как на рис. 4.
Рис. 7. Стационарное разностное решение при использовании интерполяции (4). Распределение температуры и прочие параметры приведены на рис. 6.
w
1+1/2 _ po^ l I к
i+1/2
Y — l 2
Здесь при определении второй компоненты вектора Ь используется линейная интерполяция с неопределенным пока коэффициентом а. Шаг по времени заканчивается следующим образом:
pin+1
тт
pi hx 2
po(l — к) I agpo I (l — а)gкpo
hx
(кі-1 — кі),
(pUTi — (pu)i = Po —
т lI к г-1 iN
(кг 1 — кг),
_ л>П+1 і л,П
wn+1 wn = т (Un+1/2 Un+1/2) I TPi + pi g
wi — wi = (Ui— 1 /2 — Ui+1/2 ) + т-----------------------------2-g.
Найдем из условий стационарности значения неопределенных коэффициентов а и к (к немного отличается от exp (hx/A) из-за отличия разностного решения от аналитического):
р0(1 - к) + даро + д(1 - а)кр0 = 0,
hx
Po l I к
hx 2
(l — к) I gpoК = Q.
Система имеет два решения х1>2, но нас интересует только к > 0:
pghx / 1 |
к =--------И/1 I
p
pghx
P
к
а
p
к — l pghx
(4)
2
Первые три члена разложения к в ряд совпадают с разложением в ряд ехр(Л,х/Л), также а ^ 1/2 при Л,х/Л ^ 0. Стационарный профиль скорости при определении коэффициента интерполяции по (4) приведен на рис. 5.
Случай с разностным стационарным решением при неизотермическом распределении температуры сложнее, так как схема трехточечная (в двумерном случае — девятиточечная), а зависимости величин р, е, и от компонент вектора и — нелинейные. Из-за этого тяжело физически обосновать достаточно простой разностный вид зависимости сеточных величин рг, рг так, как это было сделано в (3). Без этого же определение способа интерполяции, обеспечивающего равенство нулю потоков на целых слоях по времени, сводится или к решению сложных краевых задач, или к вычислению громоздких выражений; часто схема оказывается еще и неустойчивой. На рис. 6, 7 приведены установившиеся распределения скорости для заданной зависимости температуры от координаты при использовании старой (а = 0.5) и новой (4) интерполяций. Видно, что невязка решения при старом методе интерполяции пропорциональна шкале высот, в то время как при использовании коэффициента интерполяции (4) невязка меньше почти на порядок и зависит скорее всего уже от производной ^Л/^ж.
7. Тестирование схемы
Для тестирования предложенной схемы решено несколько задач, имеющих известное аналитическое решение. Это задача о распаде разрыва (ударной трубе), задача гидростатики (решаемая методом установления), задача о стоячей звуковой волне в прямоугольном ящике и задача о малых колебаниях среды, находящейся в равновесии в гравитационном поле при наличии теплопроводностного или объемного диссипативного механизма. Результаты тестов и более подробное описание постановки задачи даны ниже.
7.1. Распад разрыва
Задача ставится следующим образом. В начальный момент времени среда разделена абсолютно жесткой непроницаемой стенкой. В каждом полупространстве вещество однородно и покоится. В момент і = 0 стенку мгновенно убирают и надо определить дальнейшую эволюцию разрыва. Аналитически задача решается методом характеристик [18]. Приведем решение. Оно автомодельно, все величины зависят только от переменной £ = х/і. Можно выделить четыре области, где характер течения качественно меняется. В порядке от области с большим давлением к области с меньшим давлением это волна разрежения, однородное течение, контактный разрыв и за ним еще одна область однородного течения, заканчивающаяся ударной волной. Параметры среды в этих областях следующие (с*2 = у/ 7Р2/Р2, 1о = Р2/РІ):
£ Є (-то, -с^] : р = Р2, р = Р2, и = 0,
_ ! £ \2/(т-1)
( 7 - 1 £ V
£ Є £*] : р = V7, р(£) = ( 2.^2---------------------——, и = 2(^2 + £)/(7 + 1)
V 7 + 1 V7/о/
2/(7—і)
£ Є £*, v : р = /ор7, р = ( 2.°с.52 _ 1 -= ) , и = V
7 + 11/7/
£ Є [V, V] : р = рз, и = V, р = рз,
£ € [V, +то] : р = р1, и = 0, р = Р1-
Считаем, что р2 > р1, поэтому ударная волна имеет положительную скорость. Здесь переменная V обозначает скорость движения фронта ударной волны. Скорость течения V и параметры среды связаны соотношениями Гюгонио на фронте ударной волны и условиями непрерывности скорости и давления на контактных разрывах:
Р1 = (7 + !)Р1 + (7 - !)Рз рз (7 - 1)Р1 + (7 + 1)Рз ’
v = \/(Р1 - Рз)(1/Рз - 1/Р1),
Р1 + V! = Рз + (V - V)2
7 - 1 2 = 7 - 1 2 ’
и(Г) = V, Р(П = Рз-
Так как задача автомодельна, в тестах изменялся только шаг по времени. Результаты теста приведены в табл. 3. Параметры задачи следующие: р1 = р2 = 1.0, Р2 = 2.0, Р1 = 1.0, = 0.02. Аналитически найденные значения параметров решения выписаны в
первой строке таблицы; тс = /((Л.ж + ) • шах{с8 + |V|}). Скорости движения фронта
ударной волны и волны разрежения определялись как скорость движения точки, где и =
0.5v.
ТаблицаЗ
Результаты тестового расчета задачи о распаде разрыва
т/тс V Рз Рз Р2 Уг
— 1.520183 1.483218 1.26438 0.835805 -1.61383
0.868 1.51976 1.48322 1.26438 0.835756 -1.61296
0.543 1.51962 1.48322 1.26438 0.835759 -1.61297
0.216 1.5198 1.48322 1.26438 0.835767 -1.61299
7.2. Гидростатика
Решалась задача о нахождении стационарного распределения параметров среды в гравитационном поле при постоянной температуре. Задача решалась методом установления. Результаты теста приведены в табл. 4. Параметры задачи: Т = 0.4, ро = 0.3, ^ = 1.3, д = -3, Л-1 = ^д/(ЯТ) = -1.172672756. Аналитическое решение задачи следующее:
Р(х) = Ро • ехр(х/Л), р(х) = ^Т^, Л= — = ^. (5)
ДТ рд с2
Т аблица4
Результаты тестового расчета задачи гидростатики
0.04 0.02 0.01
Л-1 -1.17224 -1.17257 -1.17265
7.3. Стоячая звуковая волна
Решалась задача об определении частоты колебаний стоячих звуковых волн малой амплитуды в жестком ящике прямоугольной формы. Аналитически она выражается через скорость звука в среде, размеры области (а на Ь) и моду колебаний (п, т). Полное решение, описывающее малые возмущения в волне как функцию г и £, имеет вид
2 2 2 (т П ш = п с8 — + —
8 V а2 Ь2
Р(х, у, £) = А^р • вт(^£ + ф) • сов(птж/а) сов(ппу/Ь), пт
и(х, у, £) = А-- сов(^ + ф) • в1п(птж/а) сов(ппу/Ь),
а
пп
v(x, у, £) = А— • сов(^ + ф) • в1п(птж/а) вт(ппу/6).
Таблица5 Результаты тестового расчета задачи о стоячей волне
0.04 0.04 0.04 0.04
т 6.0Е - 3 3.0Е - 3 1.50Е - 3 0.75Е - 3
ш 17.577 17.555 17.550 17.548
7^ 2.659Е - 02 1.396Е - 02 7.008Е - 03 3.642Е - 03
0.02 0.02 0.02
т 3.0Е - 3 1.50Е - 3 0.75Е - 3
ш 18.060 18.054 18.053
7^ 3.627Е - 03 1.959Е - 03 1.060Е - 03
0.01 0.01
т 1.50Е - 3 0.75Е - 3
ш 18.230 18.228
7^ 4.545Е - 04 2.590Е - 04
Таблица6 Результаты тестового расчета задачи о стоячей волне
Лж 0.04 0.04 0.04 0.04
т 6.0Е - 3 3.0Е - 3 1.50Е - 3 0.75Е - 3
ш 12.852 12.841 12.839 12.838
7^ 1.170Е - 02 8.429Е - 03 6.542Е - 03 5.566Е - 03
0.02 0.02 0.02
т 3.0Е - 3 1.50Е - 3 0.75Е - 3
ш 12.943 12.941 12.940
7й 5.495Е - 03 5.084Е - 03 4.847Е - 03
0.01 0.01
т 1.50Е - 3 0.75Е - 3
ш 12.966 1.2966
7^ 4.716Е - 03 4.665Е - 03
Тест проводился следующим образом: задавалось возмущение скорости с заданной модой, уравнения интегрировались по времени на 10.. .20 периодов колебаний, после чего
определялась зависимость амплитуды заданной моды от времени с помощью преобразования Фурье. Полученная зависимость приближалась функцией А(і) = А0 ехр(—7^) 8Іп(ші + п). Значения 7^ и ш приведены в табл. 5 (Т = 0.4, р0 = 1.0, ^ = 1.3, 7 = 5/3, а = Ь = 1.0, п = т = 2, шан = 18.348) и табл. 6 (добавлена слабая вязкость для тестирования диссипативной части схемы: п = 2, т = 0, п = С = 1е — 4, шан = 12.9706, 7^ан = 4.60582е — 03).
7.4. Стоячая гравитационно-звуковая волна
Решалась задача об определении частоты малых колебаний среды, расположенной в гравитационном поле и ограниченной сверху и снизу жесткими стенками. Перепад плотности ртах/рт;п — 109. Для тестирования тепловой части схемы вместе с адиабатическим случаем рассмотрен случай затухающих из-за объемных потерь или теплопроводности колебаний. Частота колебаний определялась таким же образом, как и в предыдущем тесте. Аналитическое решение задачи с учетом диссипативных механизмов можно получить, линеаризовав систему уравнений (1) относительно стационарного решения (5). Система уравнений на возмущения выглядит следующим образом (для простоты изложения приводим решение для одномерной задачи, так как характерные параметры среды сильно изменяются только вдоль вектора g):
д£р + д^Рж д* + д* = ,
д^Рж д^Р
— +17 = д4р'
д ^Р 2 д ^Рж , , ^ д % , л ^ ыт
!Т +с' 1Т +<7 - 1)_дГ =<7 - ^ +(7 - ^
д£Т д / ^^р ^Р0(ж)
X дх ^дж 1чр0(ж)Д р2(ж)Я р
Здесь использовано обозначение Рж = ри и предполагается, что вектор g = еж • д. Если наложить дополнительное условие на коэффициент теплопроводности и мощность объемных потерь, то задача становится линейной (с учетом зависимости р0(ж), Ро(ж), приведенной в
(5)). Эти ограничения следующие:
х(р. Т) = рх.(Т), Ж (р, Т) = рИ-’ДТ).
Дисперсионное уравнение выглядит следующим образом:
^з + гВ^2 - (^жс2 + *7дкж)^ - гВ(к^сТ + *дкж) = 0,
^ ,2 , * ,.Л
В =(7 — 1)^ кХ + ТМ — (7 — 1)
Я \ Х А 7 'Я ^Т
При В = 0 имеем случай идеальной гидродинамики:
Т=То
, =______/Ш
Х 2А V с2 4А2'
Из условий непроницаемости стенок получаем спектр
2 2 ( 1 п2п2
Ш = ^ ( 4А2 + ~ТГ
где п = 1, 2,... — номер моды, а Ь — размер области. Сравнение численно найденной частоты колебаний с аналитическим значением приведено в табл. 7. Параметры задачи следующие: Т = 0.4, р0 = 1-0, ц = 1.3, 7 = 5/3, Ь = 4.0, п =1, ^ан = 2.02384412. Амплитуда волны за период колебаний сохраняется с точностью 10-4 % при = 0.04, 10-5 при = 0.02 и 10-6 % при = 0.01.
Теперь рассмотрим случай В = 0. Сперва находим зависимость волнового вектора от частоты и параметров системы:
; г | I 1 7^ + «7В
= -2Л у-4Л2 + 1?а ■ + «В .
Из условия непротекания на границах области снова получаем уравнение на спектр:
пп / 1 ^2 7^ + «7В
“ + в . (6)
Ь у 4Л2 с2 7^ + гВ
В случае объемных лучистых потерь В не зависит от волнового вектора и задача имеет три решения. Два из них отвечают затухающим гравитационно-акустическим волнам, а третье соответствует медленно затухающему неосциллирующему решению (которое не представляет для физики солнечной атмосферы особого интереса). Численные значения частоты колебаний и их декремента для случая объемных потерь приведены в табл. 8 (п =1, ЛЖ*/ЛТ = 1е - 2, ^ан = 2.023844040, 7ан = 2.0847514 ■ 10-4).
В случае зависимости В от волнового вектора в дисперсионное уравнение (6) надо подставить В при = —г/2Л ± пп/Ь. Но в данном случае при численном моделировании
малых возмущений возникает вопрос о постановке граничных условий для тепловой группы. При отсутствии гравитации возмущение температуры в окрестности границы было пропорционально ео8(к(ж — хгр)), что позволяло использовать условие теплоизолированной
Таблица7 Таблица8
Результаты тестового расчета задачи Результаты тестового расчета задачи о затухании
о стоячей гравитационно-акустической из-за наличия объемных потерь
волне 0.04
0.04 т 6.0Е — 3
т 6.0Е — 3 ш 2.0240234
ш 2.024003 7 2.034281Е — 4
0.04 0.02 0.04 0.02
т 3.0Е — 3 3.0Е — 3 т 3.0Е — 3 3.0Е — 3
ш 2.023987 2.023884 ш 2.0240066 2.0238888
0.04 0.02 0.01 7 2.0531185-4 2.074668Е — 4
т 1.50Е — 3 1.50Е — 3 1.50Е — 3 Лж 0.04 0.02 0.01
ш 2.023983 2.023880 2.023854 т 1.50Е — 3 1.50Е — 3 1.50Е — 3
0.02 0.01 ш 2.0240028 2.0238841 2.0238529
т 0.75Е — 3 0.75Е — 3 7 2.062304Е — 4 2.076656Е — 4 2.088439Е — 4
ш 2.023879 2.023853 0.02 0.01
т 0.75Е — 3 0.75Е — 3
ш 2.0238826 2.0238520
7 2.077350Е — 4 2.088261Е — 4
Т аблица9 Результаты тестового расчета задачи о затухании из-за наличия теплопроводности
0.04
т 6.0Е — 3
ш 2.0240146
7 1.954377Е — 4
0.04 0.02
т 3.0Е — 3 3.0Е — 3
ш 2.0239962 2.0238857
7 1.972468Е — 4 1.983270Е — 4
0.04 0.02 0.01
т 1.50Е — 3 1.50Е — 3 1.50Е — 3
ш 2.0239912 2.0238810 2.0238542
7 1.979750Е — 4 1.985893Е — 4 1.984741Е — 4
0.02 0.01
т 0.75Е — 3 0.75Е — 3
ш 2.0238796 2.0238529
7 1.987212Е — 4 1.985124Е — 4
стенки (нулевой поток тепла на границе). При наличии гравитации возмущение температуры в окрестности границы теперь пропорционально ео8(к(ж — хгр)) • ехр((х — хгр)/2Л), что делает невозможным использование граничного условия теплоизолированной стенки из-за ненулевого градиента температуры на границе. Но для малого коэффициента теплопроводности согласие должно быть более менее хорошим. Результаты такого расчета приведены в табл. 9 (п =1, х* = 1е — 2, ^ан = 2.023844046, 7ан = 2.00269867 • 10-4). Для малого выбранного коэффициента теплопроводности влияние граничного условия невелико и невязка составляет величину порядка 0.8%. Тем не менее для сравнения в табл. 10 приведены результаты аналогичного исследования при д = 0 (п =1, х* = 1е — 2,
Т аб л и ц а 10
Результаты тестового расчета задачи о затухании из-за наличия теплопроводности стоячей акустической волны
0.04
т 6.0Е — 3
ш 1.6215930
7 1.289536Е — 4
0.04 0.02
т 3.0Е — 3 3.0Е — 3
ш 1.6214995 1.6216996
7 1.281455Е — 4 1.279414Е — 4
0.04 0.02 0.01
т 1.50Е — 3 1.50Е — 3 1.50Е — 3
ш 1.6214944 1.6216945 1.6217447
7 1.276795Е — 4 1.278322Е — 4 1.281095Е — 4
0.02 0.01
т 0.75Е — 3 0.75Е — 3
ш 1.6216932 1.6217434
7 1.277666Е — 4 1.280599Е — 4
иан = 1.62175952, 7ан = 1.283429 • 10 4), когда граничные условия полностью соответствуют аналитическому решению.
8. Прогрев хромосферы опрокидывающимися акустическими волнами
Для проверки применимости описанной схемы решалась задача из области физики атмосферы Солнца — прогрев атмосферы Солнца опрокидывающимися акустическими волнами. В качестве характерных значений параметров плазмы атмосферы, используемых для обезразмеривания системы уравнений, выбраны следующие:
107
см,
Р =10"
р = 105 эрг/см3, £ = [р]/[р] = 1011 эрг/г,
V = \/Й = • 105 см/с, £ = [г]/М = л/Ш • 10 с,
д = И/М2 = 104 см/с2, п = [р][v][г] = -\/10 • 106 г/(см • с).
Т = 104 к,
В этих единицах безразмерный вид системы уравнений (1) не отличается от размерного, за исключением уравнения состояния (Л = 8.31434 • 107ерг)/(К • моль) — газовая постоянная, тильда обозначает обезразмеренные величины):
р = 10-7рЛТ/р,.
На нижней границе источник волн моделируется поршнем, скорость которого зависит от времени как и = ирвт(и£). В расчете ир = 100 м/с, и = 2п/100 с-1. Это характерные параметры источника акустических волн [3, 4]. На верхней границе стоит условие
Рис. 8. Опрокидывание генерируемой поршнем звуковой волны при распространении в стратифицированной атмосфере (в безразмерных единицах).
Рис. 9. Сформировавшаяся сильная ударная волна, распространяющаяся по начальному изотермическому фону (в безразмерных единицах).
6.0Е+0
4.0Е+0
-4.0Е+0---1—|—|—|—|—|—|—|—|—|—|—|—|—|—|—|—|—|—|—|
О.ОЕ+О 1.0Е+1 2.0Е+1 З.ОЕ+1 4.0Е+1 5.0Е+1
Рис. 10. Характерный профиль скорости для квазистационарной атмосферы, прогреваемой цугом ударных волн (в безразмерных единицах).
1.0Е+0 1.0Е-1 1.0Е-2 1.0Е-3 1.0Е-4 1.0Е-5 1.0Е-6 1.0Е-7
0.0Е+0 1.0Е+1 2.0Е+1 З.ОЕ+1 4.0Е+1 5.0Е+1
Рис. 11. Характерный профиль плотности для квазистационарной атмосферы, прогреваемой цугом ударных волн (в безразмерных единицах).
0.0Е+0
0.0Е+0 1.0Е+1 2.0Е+1 З.ОЕ+1 4.0Е+1 5.0Е+1
Рис. 12. Характерный профиль температуры для квазистационарной атмосферы, прогреваемой цугом ударных волн (в безразмерных единицах).
выпуска волн, по горизонтали стоят периодические граничные условия. Задача решалась методом установления. Начальное стационарное решение (в отсутствие источника волн) имеет следующие параметры: Т = 4000 К, |и| < 10-15, р(Х) = 0.3 • ехр(—1.07053Х).
При включении в момент времени £ = 0 источника волн от левой границы начинает распространяться звуковая волна, амплитуда которой растет по экспоненциальному закону. Потом эта волна опрокидывается и вместо звуковой волны в среде распространяется цуг ударных волн (рис. 8). Из-за низкой (по солнечным масштабам) температуры среды шкала высот мала и амплитуда ударной волны достигает значений порядка 3.68 • 103с5,
где cs = \JyRT/^ Iт—4ооок= 103 • 106см/с (рис. 9). В результате нагрева серией ударных волн шкала высот становится заметно меньше (рис. 11) и после прохождения порядка 5.. .8 ударных волн атмосфера приходит в квазистационарное состояние, в котором нагрев ударными волнами компенсирует лучистые потери. При этом скорость на фронте ударной волны после опрокидывания практически не меняется и имеет значение порядка 1.26 • 106см/с (рис. 10). Распределение параметров прогретой атмосферы, зависящих от параметров источника акустических волн, приведено на рис. 10-12.
9. Заключение
В работе обоснован выбор приближения одножидкостной газодинамики с учетом процессов лучистого теплопереноса и теплопроводности для описания течений плазмы в солнечной атмосфере. Выписана соответствующая система уравнений, построена явная консервативная схема для ее решения, исследован вопрос о погрешностях стационарного решения, полученного методом установления. Используемый алгоритм интерполяции объемных сил изменен таким образом, чтобы уменьшить погрешности при поиске стационарных решений этим методом. Применимость построенной схемы иллюстрируется кроме ряда тестов решением задачи прогрева солнечной атмосферы опрокидывающимися ударными волнами. Меньшие по сравнению с неявными схемами ограничения на диапазон изменения параметров в расчетной области плюс учет теплопроводности позволили смоделировать прогрев атмосферы Солнца ударными волнами до значительно больших высот, чем в ранее проделанных одномерных исследованиях [4, см. библиографию]. Дальнейшее увеличение исследуемого диапазона высот приводит к сильному ограничению на шаг по времени требованиями устойчивости параболической части системы (5). Для преодоления этого недостатка надо или использовать неявные методы, или в части исследуемой области воспользоваться предельным переходом х ^ то (изотермический предел).
Применение явных методов делает алгоритм удобным для распараллеливания, а использование двумерной постановки задачи позволит исследовать качественно новые задачи из области физики Солнца, такие как влияние неоднородности на прогрев солнечной атмосферы, интерференция генерируемых в фотосфере волн и целый ряд других задач [22].
Авторы выражают благодарность С. В. Алексеенко, Г. И. Дудниковой и В. А. Романову за помощь при проведении настоящей работы.
Список литературы
[1] VERNAzzA J.E., AvRETT E.H., Loeser R. Structure of the Solar chromosphere.
1. Basic computation and summary of the results // Astrophys. J. 1973. Vol. 184. P. 605-631.
[2] Ulmschneider P., Kalkofen W. ET AL. Acoustic waves in the solar atmosphere. I. The hydrodymanic code // Astronomy and Astrophys. 1977. Vol. 54, P. 61.
[3] Ulmschneider P., Kalkofen W. Acoustic waves in the solar atmosphere. III. A Theoretical tempurature minimum // Astronomy and Astrophys. 1977. Vol. 57. P. 199.
[4] Ulmschneider P., Schmitz F., Kalkofen W., Bohn H.U. Acoustic Waves in the Solar Atmosphere. V. On the chromosphere temperature rise // Astronomy and Astrophys. 1978. Vol. 70. P. 487-500.
[5] NAKAGAwA Y., Steinolfson R.S. Dynamical response of the solar corona. I. Basic formulations // Astrophys. J. 1976. Vol. 207. P. 296-299.
[6] Steinolfson R.S., Nakagawa Y. Dynamical response of the solar corona. II. Numerical simulations near the sun // Astrophys. J. 1977. Vol. 207. P. 300-307.
[7] ShibATA K., NozAwA s., MatsumOTO R. ET AL. Emergence of solar magnetic flux from the convection zone into the photosphere and chromosphere // Astrophys. J. 1990. Vol. 351. P. L25-L28.
[8] Прист Э.Р. Солнечная магнитогидродинамика: Пер. с англ. М.: Мир, 1985. 592 с.
[9] Liska R., Wendroff В. Composite schemes for conservation laws: Technical Report LA-UR 96-3589. LANL, Los-Alamos, 1996. [SIAM J. Numer. Anal. 1997].
[10] Ленг К.Р. Астрофизические формулы: руководство для физиков и астрофизиков. Ч. 1. М.: Мир, 1978. 448 с.
[11] McWhirter R.W.P., Thonemann P.C., Wilson R. The heating of the solar corona: a model based on energy balance // Astronomy and Astrophys. 1975. Vol. 40. P. 63-73.
[12] Cox D.P., Tucker W.H. Ionization equilibrium and radiative cooling of a low-density plasma // Astrophys. J. 1969. Vol. 157. P. 1157-1167.
[13] Cox D.P., Daltabuit E. Radiative cooling of low-density plasma // Astrophys. J. 1971. Vol. 167. P. 113-117.
[14] Unno W., Spiegel E.A. The Eddington approximation in the radiative heat equation // Publ. of the Astronomical Society of Japan. 1966. Vol. 18, No. 2. P. 85-95.
[15] КотЕльников И.А., Ступлков Г.В. Лекции по физике плазмы. Новосибирск: Изд-во НГУ, 1996. 128 с.
[16] Cox A.N., Stewart J.N. Rosseland opacity tables for population I compositions // Astrophys. J. Supplements Series. 1970. Vol. 19, No. 174.
[17] ЗЕльдович Я.Б., Райзер Ю.П. Физика ударных волн и высокотемпературных гидродинамических явлений. М.: Наука, 1966.
[18] Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. Т. 6. Гидродинамика. М.: Наука, 1986.
[19] Дегтярев Л.М., Фаворский А.П. Потоковый вариант метода прогонки для разностных схем с сильно меняющимися коэффициентами // ЖВМиМФ. 1969. Т. 9, №2. С. 211-218.
[20] BOUKADIDA Т., LeRoux A.-Y. A new version of the two-dimensional Lax-Friedrichs scheme // Math. Comp. 1994. Vol. 63. P. 541-553.
[21] Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989. 432 с.
[22] Alekseenko S.V., Dudnikova G.I., Romanov V.A. et al. Computational simulation of the low chromosphere heating by the shock waves’ series // Intern. Conf. on the Methods of Aerophysical Research. 2002. Pt II. P. 3-7.
Поступила в редакцию 17 января 2002 г., в переработанном виде —10 сентября 2002 г.
