Математическое моделирование ионосферы с учетом ее трехмерной неоднородности Текст научной статьи по специальности «Математика»

УДК 550.338

С. А. Ишанов, С. В. Клевцур

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ИОНОСФЕРЫ СУЧЕТОМ ЕЕ ТРЕХМЕРНОЙ НЕОДНОРОДНОСТИ

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

Difference methods an approximate solution of a three-dimensional equation of ions diffusion with mixed derivatives and divergence first derivatives are considered. The task is reduced to a chain of three two-dimensional parabolic equations on base of the summary approximation method. Difference operators are constructed, and an analysis of its properties, test calculations and an estimate of a mixed derivatives role are carried out.

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

Key words: difference method, three-dimensional equation, mixed derivatives, parabolic equation, difference operator.

Мощным инструментом исследования ионосферы Земли является вычислительный эксперимент. Основой этого эксперимента служит математическая модель ионосферы.

Задача математического моделирования структуры, динамики и теплового режима ионосферной плазмы в несогласованном виде, то есть когда информация о параметрах верхней нейтральной атмосферы задается из независимых источников (в виде некоторых эмпирических или полуэмпирических моделей), сводится к интегрированию для каждой ионизированной компоненты соответствующей системы уравнений непрерывности, движения и теплового баланса. Возьмем, например, уравнение непрерывности

^ + div(Nivi) = Qi-LiNi, (1)

at

где Ni и V — концентрация и скорость движения частиц сорта i; Qi и Li — скорость образования и коэффициент, характеризующий рекомбинационные процессы.

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

Следуя [1], построим для трехмерного уравнения диффузии аддитивные разностные схемы, обладающие суммарной аппроксимацией. Сведение трехмерной задачи к последовательности двумерных задач рассмотрим на примере уравнения диффузии (1), записанного в дивергентной форме в сферической системе координат в шаровом слое:

Щ=^(р щ ді дг I гг дг

8Ы;

г 1 1 5© ^ 00 5©

■ВД I + — (ргв — к 1 1 дг I д®

д№'

- —I Р 9© I 0Г Эг

дЫ,

— І р0і 5@1 5А.

а

5Л/;

+ аД?М5 5©

+ РвМі | +

дЫ/

+ ЗА, I?ъ ' дг

-цм, + <2,

где í — время; г — координата вдоль радиуса вектора; 0 — коширота; А — долгота; Ргг,Рхв — коэффициенты дифференциального оператора, относящиеся к иону сорта і (в данном случае

ионы О+), выражение для которых приведены в работе [2].

В дальнейшем будем опускать индекс і, относящийся к компонентам ионов. Будем искать решение Л/(£, г, ©, X) уравнения (2) в области Є (параллелепипеде), определяемую неравенствами Г0 ^ Г ^ у-уг ©О '$'■ 0 ^ 01 f Х>0 ^ А. ^ с і ра 11 її і к' її Г. ) І о 11 о.) 111 її VI урэ.внбниб (2) нэ.чэ.льными условиями (í = 0)М(0, г, ©, X) = у\і(г, ©, X) и іраничньїми условиями первого рода г, ©, ^)|г = г, ©, X) при (г, ©Д) є Г.

Функцию |а(£, г, ©, X) при г, ©Д є Г считаем непрерывной. По некоторым переменным зададим условия периодичности, например

Щї, г, ©, X) = Щ, г, ©, Х + Тх), где Т — период по долготе. Запишем (2) в операторной форме

дЫ

яы =—-гы-р = о, ді

(3)

где Ъ — дифференциальный оператор уравнения (2) по пространственным переменным; а Р = -ЬЫ + <¿2 (приЬ^0п()^0).

3

Представим оператор Z в виде суммы трех операторов Z = ^ , где

р=і

дт

Z2 =------ р

2 5©!

дт

д

дт

дт

д

2з дХ дХ+Р'к)+ дХ [Ріг дг]+ дг [РгХ дХ ^ Ргг > О, Р@0 > 0, Ри > 0.

Уравнение (2) перепишем в вице

т 1 ВЫ

= 0, %М = -—-2рМ-Рр, т = 3, (й тді

где Рр (Р ~1,т) — непрерывные функции, удовлетворяющие условию нормировки Р = ЕРр- На

Р=1

отрезке 0 ^ I ^ Т введем равномерную сетку ют = = р, ] = 0, /} с шагом т. Каждый интервал (£;> tj+1)

разобьем на т частей узлами £у+р¡т = ^ + рт /т, р = 1 ,т-1. Уравнению (2) поставим в соответствие цепочку двумерных уравнений

^ ддг ___

—— =грМ15^-Рр/ Р = 1 ,т, ^ е Ар = (^+(р_1)/т < * <

Щ(0, г, ©, %) = ф, ©, X),

^р(^+(р-1 )/т> Г’ ®' _ ^р-і(г' ®/ ^+(Р-1)/т)-

т

Краевые условия для Np задаются на части Гр, состоящей из точек пересечения Г со

всевозможными плоскостями А = const ДЛЯ Р = 1, г = const ДЛЯ Р = 2, 0 = const ДЛЯ Р = 3 и проходящими через любую внутреннюю точку (г, 0, А.) е G, ТО есть Np = \x(t, г, 0, X) при

(г, 0, X) е Гр.

За решение задачи на отрезке [f,, £,+1] принимается функция

N(tj+1, г, 0, = Nm(tj+1, г, 0, А.) = N3(tj+1, г, 0, А.).

Аппроксимируем каждое уравнение из (4) на полуинтервале t

І+( 3-і)/™

< i s= i,

-P/m ДВУХСЛОЙНОЙ

неявной разностной схемой на сетке с шагами 1гг, 1г&, 1гх и получим цепочку двумерных схем (и^/т = ЛрИ;+(3/т + ф£+|3/ш, р = 1, 2, 3, (5)

и(0, г, 0, X) - 0, ^), (г, 0, Х)&ь\, (6)

(7)

м;+Р/т =(а;+Р/т ПрИ (г,©Д)єуй/р, 7=0,/,

где узлы (г,0,/.)су/(р лежат на Гр, а Лр — разностный оператор, аппроксимирующий соответствующий дифференциальный оператор Zp. При этом и1, у1 — сеточные аналога функций N^1, г, 0, X) и Р(^, г, 0, на щ + сот.

Способы построения разностных операторов Лр рассмотрим на примере оператора Лх, аппроксимирующего Ъ\. Выражения для Л2 и Л3 могут быть получены аналогичным образом.

Разобьем сначала на составляющие = Z,т +7.г + +^г, где

^ гг - Ргг ~

дг\ дг

Z =______

г Л "г

дг

Zr& - -

dr

P

Г0

5©

д©

0r

dr

Пусть индекс г нумерует узлы по координате г, индекс к по координате 0, индекс п по координате А. Введем следующие обозначения:

р7+ Р-1 іm _ р

j+ 6-1 im ik

pj+|3/m _ p rr,i,k rr,i,k’

и так далее.

«г * = Ui,k ~ Ui-\,k Я

WMa 'К

u+\,k ui-\,k

ПК.

Аналогично обозначаются разностные отношения по другим переменным. В данном случае структура операторов Zp такова, что разностные операторы можно записать на семигочечном

шаблоне. Рассмотрим два основных шаблона, используемых при наличии смешанных производных. Показано, что оператор

\@и- [(Р^)^0 )г + (Рт&и&)? ]/2 _

Рт& ,і+1,к- 1/2иі+1, к~Рт ,0,і+1,к-1/2 иі+1,к-1~ Рт& ,і, к- 1/2иік +

+ Рт@,і,к-1/2иі,к-1+ Рт@,і,к+1 /2иі,к+1~ Ртв,і,к+1/2иік ~

1

2М©

Pi® ,i-1, k+ 1/2Ui- 1Д+ 1+ Pr& ,i-1, k+ 1/2Ui- 1Д

(8)

аппроксимирует Zr0 со вторым порядком Лг@и = Zr@M +0 |/г, |2 +|he Таким же образом записывается

A qtu — [(P0rw^ )@ + (P@rwr )q ] / 2.

(9)

Наряду с (8) можно построить оператор Лг0 на шаблоне

Лг@м = [ Сг0м0 + ^г®м0 / 2 (Ю)

со вторым порядком аппроксимации. Используя центральные разности для Лг0, получим выражение со вторым порядком аппроксимации:

Ргв,і+1,к(иі+1,к+1 иі+1,к-і) Ргв,і-1,к(иі-1,к+1 иі-1,к-і)

2 hrh@

Расчеты на контрольных примерах для модельного уравнения с известным аналитическим решением показали, что выражение (11) значительно уступает в точности выражению (8). Это вызвано прежде всего тем, что в формуле (11) не захватывается центральный узел (і, к). В связи с этим дальнейший анализ производится только для шаблона 8.

Операторы Лгг и Л г записываются известным образом [3]. В результате разностный двухмерный эллиптический оператор получит вид

ВцЦі.к-1 + ЦїА'+ІД-І + ^гіД'-ІД " СцЦік + Ц'іД'+ІД + ЦіДг-ІД+І + ^яДг'Д+І / (л ^

і = ТЩ, к = ЇЖк.

Одним из наиболее эффективных способов решения системы разностных уравнений (12) является итерационный а —(3 алгоритм [3], который представляет собой двумерный вариант метода прогонки по двум пространственным координатам. Условием его сходимости является аналог диагонального преобладания по строкам у матрицы системы разностных уравнений в одномерном случае, записываемом в виде

Сік 5= Вік + ^ік + Кік + Еік + Вік + Уік. (13)

Условие (13) сводится к неравенству

~[Рг,і+1/2,к ~Рг,і-1/2,к]/К +(1ік +1/т0, (14)

где величина сіік является сеточным аналогом функции Ь(1^р/тгіг®к,%п).

В ионосферных приложениях неравенство (14) требует для своего выполнения практически неприемлемых малых шагов т по времени. Из (14) также следует, что наличие смешанных производных не влияет на «диагональное» преобладание и качество разностной схемы определяется наличием в операторе составляющей 2Г — первой производной дивергентного вида. В случае трехмерного уравнения (2) возникают дополнительные аналогичные осложнения, связанные с добавлением составляющих и причем о пространственно-временном поведении коэффициентов Р0((, г, ©Д) и Р,(1, г, 0, /.) ничего не известно. В связи с этим возникает проблема построения безусловно устойчивых разностных схем при наличии производных дивергентного вида по переменным г, ©, А,. Следуя [3], рассмотрим преобразование оператора

€гг + ^ = — ^п-

гг Р

цг = ехр \y-dr'. (16)

г0 п

Разностные аналоги для (15) строятся с помощью интегро-интерполяционного метода [1] на шаблоне (8). В результате получаем разностный оператор (12), но с «диагональным» преобладанием по столбцам:

Сік > Ві к+1 +Ьік + Кі+1к +ЕІ_1 д + Оік + Уік_г. (17)

Для этого случая нами разработан модифицированный итерационный а —Р алгоритм, сходимость которого обеспечена безусловным (в отличие от неравенства (14)) выполнением неравенства (17). Этот алгоритм протестирован на ряде контрольных примеров уравнений, имеющих аналитическое решение. На основе ряда трехмерных математических моделей ионосферной плазмы в шаровом слое, в состав которых входил данный алгоритм, рассчитаны реальные геофизические ситуации [3]. В частности, проведена оценка роли смешанных производных в уравнениях диффузии. Показано, что их неучет может привести к 15-процентному отклонению в электронной концентрации Б2-слоя ионосферы на средних широтах в околополуночные часы.

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

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

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

—+ргм! = —

дг ) дт

Р д

±ГГ ^

Уг &

ігМ]

(15)

1. Самарский А. А. Теория разностных схем. М., 1977.

2. Фаткуллин М. Н., Клевцур С. В., Латышев К. С. Оператор переноса в уравнениях непрерывности для ионов в трехмерно неоднородной области F // Геомагнетизм и аэрономия. 1984. Т. 24, № 6. С. 906 — 910.

3. Игианов С. А., Клевцур С. В., Латышев К. С. Алгоритм а — Р итераций в задачах моделирования ионосферной плазмы / / Математическое моделирование. 2009. 21, № 1. С. 33 — 45.

4. Ишанов С. А., Клевцур С. В., Латышев К. С. и др. Многопроцессорная реализация одного итерационного алгоритма для систем с распределенной памятью // Вестник Российского государственного университета им. И. Канта. 2009. Вып. 10. С. 64 — 73.

Об авторах

С. А. Ишанов — канд. физ.-мат. наук, доц., РГУ им. И. Канта, e-mail: math@dekan.albertina.ru

С. В. Клевцур — канд. физ.-мат. наук, доц., РГУ им. И. Канта.

Authors

Dr. C. A. Ishanov — assistant professor, IKSUR, e-mail: math@dekan.

albertina.ru

Dr. C. V. Klevtsur — assistant professor, IKSUR.