Модель среднемасштабных неоднородностей низкоширотной ионосферы Земли: оценка метода нелинейной коррекции Z-схемы Текст научной статьи по специальности «Математика»

Об авторе

Юрий Иванович Попов — канд. физ.-мат. наук, проф., Балтийский федеральный университет им. И. Канта, Калининград. E-mail: [email protected]

About the author

Dr Juriy Popov — Prof., I. Kant Baltic Federal University, Kaliningrad. E-mail: [email protected]

УДК 550.388.2

Н. М. Кащенко, С. В. Мациевский, А. А. Викторов МОДЕЛЬ среднемасштабных неоднородностей

НИЗКОШИРОТНОЙ ИОНОСФЕРЫ ЗЕМЛИ: ОЦЕНКА МЕТОДА НЕЛИНЕЙНОЙ КОРРЕКЦИИ Z-СХЕМЫ

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

A nonlinear correction of the finite-differential scheme for solution of the cross convection-diffusion equations within models of the Rayleigh — Taylor instability in the equatorial region of the Earth ionosphere is considered. For test problems the experimental value of the approximation order of the offered method of non-linear correction of the finite-differential scheme for different types of limiters is received numerically.

Ключевые слова: математическое моделирование, численное моделирование, уравнение переноса, монотонная разностная схема, ионосфера.

Key words: mathematical modeling, numerical simulating, convection-diffusion equations, monotonic finite-differential scheme, ionosphere.

11

1. Математическая модель экваториального Б-слоя ионосферы Земли

Хорошо известно, что Рэлей-Тейлоровская неустойчивость в экваториальной ионосфере Земли приводит к развитию плазменных неоднородностей с характерными масштабами по времени 10—1000 с и по пространству 1 — 100 км [1; 4]. При математическом исследовании таких процессов могут использоваться в соответствии с работами [1 — 8] следующие два приближения:

1) плазма тепловая;

2) плазма равновесная.

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

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

12

Поэтому будем рассматривать ионосферу Земли как сплошную среду и воспользуемся для ее моделирования уравнениями магнитной гидродинамики. Пусть также ионосферная плазма Б-области обладает упрощающими свойствами:

1) плазма замагничена;

2) плазма квазинейтральна;

3) выполнены условия электростатики, и следовательно, потенциальности электрического поля:

Е = -УЫ, (1)

где и — электрический потенциал; Е — напряженность электрического поля;

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

—и к 2 1

ди.

дТ

( дТ. ^ ^ ^

+ р)У1У) +У±Ц , = С. - Р; , (3)

\

+ у(иУ) = Q -Ц, (2)

(Уу1 )Т

да 11

/

У1 (< У1и) = У1 А, (4)

где и, У1, Qj, Ц, р.Т., ц ., С., Р. — соответственно концентрация, дрейфовая скорость, скорости образования и потерь заряженных частиц, давление, температура, плотность теплового потока, скорость нагрева и скорость охлаждения частиц сорта к — постоянная Больцмана; < — тензор интегральной проводимости; А — интегральные источники тока. Характерной особенностью модели (1) —(4) является выполнение условия несжимаемости сплошной среды:

У1(у.) = 0. (5)

Полученные путем таких упрощений двумерные уравнения решаются в экваториальной плоскости Земли. Для вычисления коэффициентов уравнений системы использована глобальная термосферная модель МБК [9].

Кроме того, для решения этих нестационарных уравнений необходимо задать начальные значения, содержащие при необходимости начальные возмущения.

Зададим простую область решения в виде прямоугольника в координатной системе (у, г), где у — горизонтальная координата, а 2 — вертикальная.

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

Ц- ц =

(6)

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

Наконец на боковых границах прямоугольной области интегрирования заданы условия

ди■ дТ

= 0,-!- = 0. (7)

дУ дУ

2. Численная модель экваториального Г-слоя ионосферы Земли

Описанная в п. 1 математическая модель (1) — (4) представляет систему, состоящую из:

1) уравнений поперечного к магнитному полю Земли переноса заряженных частиц и температуры;

2) уравнения потенциала электрического поля, имеющего эллиптический тип.

Теперь займемся исследованием возможности применения в рассмотренных задачах метода решения уравнений переноса, построенного на основе проанализированной в работе [10] разностной схемы для решения одномерного уравнения конвективного переноса. Данная схема (названная в [10] «^-схемой») имеет второй порядок аппроксимации по времени и координате и абсолютно устойчива. Однако в силу теоремы Годунова эта схема немонотонна.

Для тестирования этой схемы рассмотрим двумерное уравнение переноса, записанное в следующем виде:

ди т, ди т, ди — + К — + Уг—

дt у ду дг

+ Уу^ + У^ = 0, (8)

где (Уу, Уг) — заданное поле скоростей, для которого выполнено условие несжимаемости; и — одна из неизвестных (щ или Т).

Поскольку операторы переноса по у и г некоммутативны, то для решения такой двумерной задачи использовалась аддитивно симметрированная схема расщепления, которая для решения двумерного уравнения (8) на отрезке времени [Ь0, Ь0 + т], где т — шаг по времени Ь, имеет вид

и + уу Ёи

дЬ у ду

+ Уу = 0, и,(М = и(М, Ь 6 [*0, *0 +т],

^ + у ^ = 0, í 6 [^ +т, Ь0 + 2т],

дЬ дг

^ + Уг ^ = 0, иа(М = и(М, Ь 6 [*0,10 +т], (9)

дЬ дг

^ + уу Ё^ = 0, * 6 [Ь0 +т, Ь0 + 2т],

дí у ду 0 0

1

и(^0 + т) = (*0 +2т) + иъ (^ +2т)).

13

14

Как известно, такая схема имеет точность по времени 0(т2), если каждый шаг схемы имеет точность 0(т2).

Далее, очевидно, что каждый из шагов нашей схемы расщепления (9) можно записать в виде

^ + V ^ = 0, dt dx

(10)

где V — скорость переноса; t — время; x — координата.

Теперь будем считать сетку равномерной, а шаги сетки по переменным t и x обозначим соответственно через т и h. На рисунке 1 показан шаблон этой разностной схемы для условий V = const и V > 0, используемой в уже упомянутой работе [10].

Рис. 1. Шаблон разностной схемы для одномерного уравнения

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

Исходная разностная схема согласно [10] имеет вид

nf+1 - nf 1

f „„m+1 „„m+1

n

- n-1.+—- n-

h

h

= 0,

(11)

где верхние индексы — это номера узлов по времени, а нижние — номера узлов по пространственной переменной.

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

,„m+1 „„m+1 „„m+1

П; — П; П; — П; i

+ Vj-hJ—L +

nm+1 - nf+1 ——-— fi 1/2V 2h J--1/2

nf - nf+11 2h

= 0, (12)

где / с индексами — это корректирующие множители.

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

f -+1/2 = f

f „„f -.„f+1 Л n- - ni-1 -„f -„f+1 v nm+1- nm у

= f (r),

(13)

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

х

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

Целью данной статьи стало описание результатов численных экспериментов, проведенных для сравнения разных типов ограничителей, монотонизирующих разностную схему (12).

Для приводимых ниже экспериментов выбраны следующие четыре ограничители, взятые из работ [11 — 14]:

1) ограничитель № 1 "МтМод":

/(г) = тт(г, 1), г > 0, /(г) = 0, г < 0;

2) ограничитель № 2 "БирегВее":

\/(г) = тах(тт(2г, 1), тт(г, 2)), г > 0, [/(г) = 0, г < 0;

3) ограничитель № 3 Ван Лира:

\/ (г) =

2г

, г > 0,

4) ограничитель № 4:

1 + г [/(г) = 0, г < 0;

а \ гл/2

/(г ) = I „■, г > 0,

41-

+ г

/(г) = 0, г < 0.

На рисунке 2 приведены графики этих четырех ограничителей. Как известно [8 — 11], область монотонности ограничена графиками МтМод и БирегВее.

15

Рис. 2. Графики ограничителей, использованных в численных экспериментах

16

3. Результаты численных экспериментов

Область интегрирования для уравнения (8) задавалась в виде

-500 км ^ у ^ 500 км, 100 км ^ z ^ 1100 км. Профиль скоростей задавался модельно в виде

Уу = 1,05Ут + Ут(2 - 2()), У2 = 1,35Ут - Ут(У - У0), (14)

" г0 г0

где Ут = 100 м/ с, г0 = 500 км.

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

(ус - 339, 2с - 350), (ус - 350, 2с - 450), (ус - 250, 2с - 460).

При этом основание высоты профиля соответствует точке

(ус - 313, 2с - 420).

В этих формулах точка (ус, 2с) — центр области решения.

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

Р =

1(п - и0)2

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

Для решения с описанным начальным профилем и заданным полем скоростей проведено экспериментальное определение порядка аппроксимации. Для этого были выполнены расчеты с разными шагами по пространству и соответственно по времени, результаты которых для момента времени 4160 секунд от начала решения показаны в таблице. Число Куранта в этих расчетах равно 0,667.

Относительные погрешности решения с разными ограничителями

Число ячеек разностной сетки по у и г и число шагов по времени Номер ограничителя

1 2 3 4

250; 735; 550 0,131 0,0489 0,0794 0,0890

500; 750; 1100 0,0709 0,0258 0,0393 0,0444

1000; 1500; 2200 0,0371 0,0144 0,0192 0,0217

2000; 3000; 4400 0,0191 0,00835 0,00937 0,0106

4000; 6000; 8800 0,00974 0,00528 0,00457 0,00511

17

Среднее значение отношения погрешностей для этой серии расчетов соответствует следующим погрешностям:

1) 0(йа94) для ограничителя МтМо^

2) 0(йа8°) для ограничителя БирегВее;

3) 0(Й1,03) для ограничителя Ван-Лира;

4) 0(й1,03) для ограничителя четвертого вида.

При этом лучшая точность получена при использовании ограничителя Ван-Лира и немного хуже при использовании ограничителя последнего вида.

На рисунке 3 приведены профили величины и для разных моментов времени для сетки с числом узлов 251 по у и 376 по г для ограничителя типа Ван-Лира. -

100 .400 О 400 V, км -400 0 400 У, ™ -400 0 400 у. ™

А) Б) В)

Рис. 3. Распределение величины п в моменты времени, соответствующие шагам с номерами 1 (А), 280 (Б) и 550 (В)

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

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

18

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

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. Ботова М. Г., Романовская Ю. В., Намгаладзе А. А. Вариации ионосферы: сопоставление результатов моделирования с данными наблюдений // Вестник МГТУ. 2014. Т. 17, № 2. С. 385-393.

7. Huba J. D., Joyce G., Krall J. Three-dimensional modeling of equatorial spread F / Aeronomy of the Earth's atmosphere and ionosphere. I AG A Special Sopron Book Series. 2011. Vol. 2. P. 211-218.

8. Krall J., Huba J. D. SAMI3 simulation of plasmasphere refilling // Geophys. res. lett. 2013. Vol. 40. P. 2484-488.

9. 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 / / Ibid. 1977. Vol. 82, № A1. P. 2139-2147.

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

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

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

13. 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.

14. Harten A. High resolution schemes for hyperbolic conservation laws. // J. of Comput. Physics. 1983. Vol. 49. P. 347-393.

Об авторах

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

E-mail: [email protected]

Сергей Валентинович Мациевский - канд. физ.-мат. наук, доц., Балтийский федеральный университет им. И. Канта, Калининград.

E-mail: [email protected]

Андрей Александрович Викторов - асп., Балтийский федеральный университет им. И. Канта, Калининград.

E-mail: [email protected]

About the authors

Dr Nikolay Kashchenko — Ass. Prof., I. Kant Baltic Federal University, Kalinin-

E-mail: [email protected]

Dr Sergey Matsievsky — Ass. Prof., I. Kant Baltic Federal University, Kaliningrad. E-mail: [email protected]

Andrey Viktorov — PhD student, I. Kant Baltic Federal University, Kaliningrad. E-mail: [email protected]

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

For the numerical solution of the boundary value problem of convective heat transfer in a ladle while blowing steel argon proposed finite-difference approximation of the original mathematical model of convective heat transfer in a continuous divergence form. Used stable and economical implicit monotone conservative difference scheme.

Ключевые слова: конвективный теплоперенос, уравнение переноса вихря, уравнение переноса тепла, уравнение Пуассона для функции тока, дивергентное представление модели, разностная схема, односторонняя четырехточечная разностная аппроксимация внутрь области.

Key words: convective heat transfer, vortex transport equation, heat transfer equation, Poisson equation for the stream function, divergent representation of a model, differences scheme, one-sided four-point difference approximation inside the area.

Математическая модель конвективного теплопереноса в ковше со сталью в дивергентном представлении может быть представлена в вице (©, ю, Т)-системы [1; 2]

grad.

19

УДК 669.18.046.517

С. В. Веревкин, С. А. Дёмин

КОНЕЧНО-РАЗНОСТНАЯ АППРОКСИМАЦИЯ КРАЕВОЙ ЗАДАЧИ КОНВЕКТИВНОГО ТЕПЛООБМЕНА В КОВШЕ СО СТАЛЬЮ

(1)

© Веревкин С. В., Дёмин С. А., 2016

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