Использование Z-схемы в модели неустойчивости Рэлея - Тейлора Текст научной статьи по специальности «Физика»

УДК 550.388.2

Н. М. Кащенко, С. В. Мациевский

ИСПОЛЬЗОВАНИЕ Z-СХЕМЫ В МОДЕЛИ НЕУСТОЙЧИВОСТИ РЭЛЕЯ - ТЕЙЛОРА

Рассмотрена нелинейная разностная схема решения уравнений переноса в рамках моделей неустойчивости Рэлея — Тейлора в экваториальной области ионосферы Земли. Для тестовых задач численно подтверждена монотонность построенной схемы. Получено экспериментальное значение порядка аппроксимации предлагаемого метода нелинейной коррекции разностной схемы.

A nonlinear finite-differential scheme for solution of convection-diffusion equation in the field of models of Rayleigh-Taylor instability in the equatorial region of Earth ionosphere is considered. For test tasks monotony of the constructed scheme is in number confirmed. Experimental value of an order of approximation of the offered method of nonlinear correction of the finite-differential scheme is received.

Ключевые слова: конечно-разностная схема, нелинейность, уравнение переноса, монотонность, неустойчивость Рэлея — Тейлора, ионосфера Земли.

Key words: finite-differentia! scheme, nonlinearity, convection-diffusion equation, monotony, Rayleigh-Taylor instability, Earth ionosphere.

© Кащенко Н. М., Мациевский С. В., 2016

Вестник Балтийского федерального университета им. И. Канта. 2016. Сер. : Физико-математические и технические науки. № 1. С. 18 — 25.

При исследовании неустойчивости Рэлея — Тейлора (НРТ) в экваториальной ионосфере используется, в соответствии с работами [1 — 5], подход с приближениями тепловой, равновесной плазмы, поэтому будем рассматривать ионосферу как сплошную среду и воспользуемся для ее моделирования уравнениями магнитной гидродинамики. Согласно Б. Н. Гершману [2], систему уравнений Максвелла и гидродинамических уравнений с учетом электромагнитных сил можно записать в виде, включающем пять уравнений: непрерывности, движения ионов и электронов, теплопроводности ионов и электронов, потенциальности электрического поля и непрерывности электрического тока [1; 4]:

^+ v ,) = Q,-L,,

~я+ w=- nm+m (ё+v,x в)-ь>v- v )* (V>-v )+g,

, , , trj

3 JdTj

—nk 2 ,

\

+(vV± )t + р, vv, +v± q, = с,- p,

v Ui /

Ух Е = 0, (1)

V = в1п1У1 =0,

где п , V,, Q;> Ь , т, е, р , у;П, V ^ Ту, у, G , Р , — соответственно концентрация, дрейфовая скорость, скорости образования и потерь, масса, заряд, давление, частоты соударений с нейтралами, частоты соударений между заряженными частицами, температура, плотность теплового потока, скорость нагрева и скорость охлаждения частиц сорта у к — постоянная Больцмана; у — плотность тока; Е — напряженность электрического поля.

Вследствие замагниченности ионосферной плазмы Б-области при использовании приближений нейтральности ионосферной плазмы, выполнении условия электростатики (1), а следовательно, потенциальности электрического поля

Е = -уи,

где и — электрический потенциал, и в приближении сильной вытяну-тости плазменных неоднородностей вдоль магнитного поля Земли можно использовать двумерное приближение процесса [1; 4], уравнения которого имеют вид

дп

И

L+v± (nV) = Q-L, (2)

-nk 2 1

& + vv )T J+p, vv +v q, =G, - P, (3)

v± (< v^) = v± A,

19

где < — тензор интегральной проводимости; A — источники тока.

20

Для вычисления коэффициентов уравнений системы могут быть использованы глобальные термосферные модели, например MSIS [6].

Для решения этих нестационарных уравнений необходимо задать начальные значения, содержащие при необходимости начальные возмущения. Уравнения (2), (3) этого приближения описывают поперечный силовым магнитным линиям перенос частиц и энергии соответственно.

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

Исследуем возможность применения метода решения уравнений переноса, построенного на основе предложенной в [7] новой разностной «Z-схемы», которая может быть использована для решения простейшего одномерного уравнения конвективного переноса.

Эта схема имеет второй порядок аппроксимации по времени и координате и абсолютно устойчива. При этом в [7] утверждается, что она монотонна, что неверно в силу теоремы Годунова.

Для тестирования рассмотрим простейшее одномерное уравнение переноса вида

^ + У ЁП = 0,

dt dx

где n — концентрация; V — скорость переноса; t — время; x — координата. При этом будем считать, что V = const и V > 0.

Обозначим шаги равномерной сетки по переменным t и x через т и h соответственно.

На рисунке 1 показан шаблон разностной схемы (Z-схемы), предложенной в работе [7].

Рис. 1. Шаблон разностной схемы (Z-схемы)

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

Рп =

Рй =

Х(п ~ по)2

Ло)2

Х л

(4)

(5)

где рп — погрешность концентрации; рй — погрешность разностной производной; п0 — точное решение; й — разности первого порядка от величины п; й0 — разности первого порядка от точного решения п0.

Суммирование проводится по всем узлам сетки для формулы (4) и по всем ячейкам сетки для формулы (5).

Для получения монотонной схемы используем подход, описанный в работах [8 — 10]. Но в отличие от этих работ для коррекции схемы используем не потоки, а их аналоги («косые потоки»).

Исходная разностная 2-схема, согласно [7], имеет вид

пт+1 ~ пт

/„„т+1 „„т+1 „„т „„т п, ~ п,~1 п,+1 ~п,

к

к

= 0.

21

Эту схему перепишем в следующем виде, одновременно добавив корректирующие множители:

пт+1 ~ пт -+ У

„„т+1 „„т+1 „„т „„т-

п. ~ п. п,+1 ~ п.

-Л_+ / У._'-

^ +1/2 у

к

2к

/.~1/2У

„„т „„т+1 п, ~ п,~1

2к

= 0,

где / с индексами — корректирующие множители. Если / = 1, то получаем схему из работы [7]. Если величины / выбрать в виде функций

//+1/2 1

( пт ~ п^ Л пт ~ пт+1

\ГЧ+1 ГЧ

= / (г),

то при надлежащем выборе вида функций /(г) получим монотонную схему.

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

/ (г) =

2г

, г > 0,

1 + г /(г) = 0, г ^ 0.

Для численных экспериментов было выбрано У = 1, а область решения х е [0, 200], t е [0, 80]. Тестирование проведено для начальных профилей, показанных на рисунке 2.

2

0

т

X

22

с й

Рис. 2. Начальные профили для п:

а — трапеция; Ь — треугольник; с— двойной треугольник; й — гладкий колоколообразный профиль

На рисунке 3 показаны профили, полученные после 1200 шагов по времени для начального профиля типа й.

Рис. 3. Профили п через 1200 шагов по времени. Справа показана увеличенная в 10 раз верхняя часть левого рисунка: тонкая линия — точное решение, толстая линия — численное решение

Ь

На рисунке 4 показана зависимость погрешности рп для профиля типа трапеции а от номера N шага по времени.

3.5 4.5 5.5 6.5 1пЫ

-5 \

-5.5

-б-

-6.5-1прп

Рис. 4. Зависимость погрешности рп для профиля типа а от номера N шага по времени

Эта зависимость хорошо аппроксимируется линейной формулой 1п рп = 0,371-1п N ~ 7,565,

что соответствует

рп = 5,18-10-4Ы0,371. Аналогичная зависимость для профиля типа Ь имеет вид 1п Рп = 0,405-1п N ~ 7,432,

что соответствует

Рп = 5,92-10-4Ы

,-4^0,405

Для задач неустойчивости важным параметром является также погрешность разностной производной. На рисунке 5 показана зависимость погрешности самой величины и ее разностной производной для гладкого начального профиля типа й от количества шагов по времени.

3.5 4.5 5.5 6.5

-4

1прп,1п(1п

Рис. 5. Зависимость погрешностей для профиля типа й от номера шага по времени: а — йп, Ь — рп

Для гладкого решения с начальным профилем типа ^ который был задан формулой

п = а((х ~ 38)(х ~ 62))2,

причем параметр а в этой формуле задан так, чтобы максимальное значение п было равно 100, проведено экспериментальное определение

23

порядка аппроксимации. Для этого были выполнены расчеты с разными шагами по пространству и, соответственно, по времени. Их результаты для момента времени, равного 80, представлены в таблице.

Экспериментальное определение порядка аппроксимации

h 0,1 0,05 0,025 0,0125 0,00625

pn 3,05-10-3 8,66-10-4 2,45-10-4 6,99-10-5 2,00-10-5

24

Обработка результатов, показанных в таблице, дает следующую оценку погрешности для этой серии расчетов: рп и 0,195-й1,81.

Приведенные результаты были получены для шага по времени, равного 2/3 от шага Куранта. При превышении шага Куранта построенная схема теряет монотонность, но остается устойчивой. Таким образом, построенная разностная схема при шагах по времени, не превышающих шаг Куранта, является монотонной и показывает хорошие точностные характеристики для решений разных типов.

Работа выполнена при финансовой поддержке РФФИ по проекту 14-01-00020.

Список литературы

1. Гайдуков В. Ю., Кащенко Н. М., Мациевский С. В. и др. Запуск экваториальных пузырей путем модификации E-слоя // Геомагнетизм и аэрономия. 1991. Т. 31, № 6. С. 1042-1048.

2. Гершман Б. Н. Динамика ионосферной плазмы. М., 1974.

3. Грэд Г. О кинетической теории разреженных газов // Механика. 1952. Вып. 4. С. 71-79. Вып. 5. С. 61-96.

4. Мациевский С. В., Кащенко Н. М., Никитин М. А. Ионосферные пузыри: ионный состав, скорости движения плазмы и структура // Известия вузов. Радиофизика. 1989. Т. 32, № 11. С. 1320-1326.

5. Рыбин В. В., Поляков В. М. Об амбиполярности движений ионосферной плазмы // Ионосферные исследования. 1983. № 33. С. 5—44.

6. Hedin A. E., Salah J. E., Evans J. V. et al. A global thermospheric model based on mass spectrometer and incoherent scatter data MSIS 1. N2 density and temperature / / J. Geophys. Res. 1977. Vol. 82, № A1. P. 2139-2147.

7. Кольцова Э. М., Федосова Н. А., Балашкина Ю. А. Новый метод разностной аппроксимации решения для задач механики сплошных сред // Успехи в химии и химической технологии. 2014. Т. 28, № 1. С. 64-66.

8. Ладонкина М. Е., Неклюдова О. А., Тишкин В. Ф., Чеванин В. С. Об одном варианте существенно неосциллирующих разностных схем высокого порядка точности для систем законов сохранения // Математическое моделирование. 2009. Т. 21, № 11. С. 19-32.

9. Сафронов А. В. Оценка точности и сравнительный анализ разностных схем сквозного счета повышенного порядка // Вычислительные методы и программирование. 2010. Т. 11. С. 137-143.

10. Van Leer B. Upwind and high-resolution methods for compressible flow: from donor cell to residual-distribution schemes // Commun. Comput. Phys. 2006. Vol. 1, № 2. P. 192-206.

Об авторах

Николай Михайлович Кащенко — канд. физ.-мат. наук, доц., Балтийский федеральный университет им. И. Канта, Калининград. E-mail: kaschtschenko@mail.ru

Сергей Валентинович Мациевский — канд. физ.-мат. наук, доц., Балтийский федеральный университет им. И. Канта, Калининград. E-mail: matsievsky@newmail.ru

25

About the authors

Dr Nikolay Kashchenko, ass. prof., I. Kant Baltic Federal University, Kaliningrad. E-mail: kaschtschenko@mail.ru

Dr Sergey Matsievsky, ass. prof., I. Kant Baltic Federal University, Kaliningrad. E-mail: matsievsky@newmail.ru