CC BY
37
Вычислительные технологии
Том 15, № 5, 2010
Расчет глобальных электрических полей в земной атмосфере*
В. В. Денисенко, Е. В. Помозов Институт вычислительного моделирования СО РАН, Красноярск, Россия
e-mail: denisen@icm.krasn.ru
Построен многосеточный алгоритм для решения глобальной задачи электропроводности атмосферы. Решена задача о проникновении электрического поля из ионосферы до поверхности Земли в спокойных геомагнитных условиях. Показано, что электрическое поле вблизи поверхности Земли можно строить как решение одномерной по высоте задачи, но чтобы получить поля и токи в верхней атмосфере, необходимо решать трехмерную задачу.
Ключевые слова: атмосфера, электрическое поле, эллиптические уравнения, параллельные вычисления.
Введение
При изучении атмосферных электрических полей основное внимание обычно уделяется грозовым явлениям. Глобальным результатом гроз является поддержание разности потенциалов между поверхностью Земли и ионосферой около 300 кВ, Поскольку атмосфера является проводником, эта разность потенциалов обеспечивает примерно постоянный электрический ток из ионосферы на землю. Соответствующее электрическое поле вблизи поверхности Земли в среднем составляет 130 В/м. Вместе с тем на фоне средней картины происходят изменения глобального электрического поля как за счет генераторов, расположенных выше атмосферы, так и за счет подземных источников тока.
Логично возникает вопрос об использовании электрических полей подземного происхождения для анализа происходящих под землей процессов и, в частности, как предвестников землетрясений. Для создания глобальной системы мониторинга актуальна перспектива измерения этих процессов с помощью космических аппаратов. При этом необходимы знания процессов проникновения полей от земли через атмосферу до ионосферы. Данному вопросу посвящены многочисленные исследования, обзор которых дается, например, в работах [1, 2]. Некоторые новые результаты в данном направлении приведены в [3].
Интересен и обратный эффект — проникновение ионосферных электрических полей через атмосферу до земли.
Для моделирования этих явлений необходима математическая модель атмосферного глобального проводника, созданию которой и посвящена настоящая работа.
*Работа выполнена при финансовой поддержке РФФИ (гранты № 07-05-00135, 09-06-91000) и программы РАН № 16.3.
Уравнение электропроводности атмосферы является эллиптическим уравнением со скалярным коэффициентом, значение которого изменяется с высотой в миллионы раз. Необходимо решать различные смешанные краевые задачи в трехмерных областях с горизонтальными размерами, на два-три порядка превышающими вертикальный размер,
1. Основные уравнения
Стационарные электрические ноля удовлетворяют законам Фарадея, сохранения заряда и Ома:
Ух Е = 0, (1)
Уj = 0, (2)
j = аЕ, (3)
где Е — напряженность электрического поля, — плотность тока, а — проводимость, В атмосфере ниже 70 км проводимость изотропна, поэтому а — скаляр. Уравнение (1) позволяет ввести электрический потенциал V такой, что
Е = -УК (4)
Тогда система (1)-(3) сводится к уравнению
-У(аУ V) = 0, (5)
а
2. Проводимость и граничные условия
Типичное высотное распределение проводимости в атмосфере Земли приведено в |4| и показано на рис, 1, От поверхности Земли вверх до 85 км проводимость возрастает на восемь порядков.
Рис. 1. Высотное распределение проводимости в атмосфере [4|; справа показано разбиение но высоте на 16 слоев
Выше, в ионосфере, проводимость становится тензором и его диагональные компоненты увеличиваются еще на несколько порядков, В основном проводящем слое ионосферы проводимость в направлении вектора магнитной индукции В в сотни и тысячи раз превышает остальные компоненты тензора проводимости. Поэтому при моделировании крупномасштабных ионосферных электрических полей принято использовать двумерные модели, в которых магнитные силовые линии являются эквипотенциальными, Такого рода модель описана в работе [5].
Оба указанных упрощения не справедливы в слое, который является переходным между атмосферой и ионосферой и обычно находится на высотах 70-90 км. Здесь с ростом высоты проводимость перестает быть изотропной, как в атмосфере, но в направ-В
как в основных слоях ионосферы. Если исключить из рассмотрения экваториальную область, где геомагнитное поле почти горизонтально, каждая магнитная силовая линия при прохождении рассматриваемого слоя, имеющего высоту около 20 км, смещается по горизонтали не более чем на несколько десятков километров. Поэтому, если экстраполировать эквипотенциальноеть магнитных силовых линий в этот слой, распределение ионосферного электрического потенциала, полученное на высоте 90 км, при переходе на высоту 70 км перераспределится по горизонтали на эти десятки километров. Такое смещение можно рассматривать как оценку погрешности, вносимой тем, что в нашей модели нижняя ионосфера заменена поверхностью, на которой проводимость скачком меняется от изотропной атмосферной до существенно анизотропной ионосферной. Естественно, данная модель применяется только для полей с горизонтальным масштабом не менее сотни километров. Проведенные тестовые расчеты показали, что выбор такой поверхности на высотах 80-90 км незначительно меняет получающееся атмосферное электрическое поле. Корректный учет этого переходного слоя может быть сделан в рамках трехмерной модели электропроводности с гиротропным тензором проводимости. Соответствующая модель проводимости приведена в [5], энергетический метод решения такой задачи, имеющей несамосопряженный оператор, предложен в [6], В этой, более адекватной, модели целесообразно отделить атмосферу, добавив условия сшивки на новой границе, В итоге для атмосферы получится модель, применяемая в настоящей работе, В силу изотропии проводимости методы решения для собственно атмосферы существенно упрощаются,
В работе также используется обычное предположение о малости влияния атмосферной проводимости на распределения электрических полей в ионосфере, сделанные в силу малости атмосферной проводимости по сравнению с ионосферной. Поэтому в предлагаемой модели электропроводности атмосферы на ее верхней границе полагается заданным распределение электрического потенциала V:
где г, в, р — сферические координаты, Яир — радиус сферы, которая используется в качестве верхней границы атмосферы.
Нижней границей атмосферы является поверхность Земли, Поскольку проводимости морской воды и влажной почвы намного превышают проводимость приземного воздуха, Землю полагаем эквипотенциальной:
V \Р== = Шр)
(6)
V\г=ЯЕ = 0
(7)
где Ке — радиус Земли, Это условие записано в предположении сферичности Земли, Если учитывать рельеф с заданной с помощью функции Н(9, ф) высотой поверхности над уровнем моря, то условие Дирихле (7) принимает вид
К 1=яЕ = (8)
Если рассматривать плохопроводящие участки суши, такие как скальные массивы, то на соответствующих частях поверхности Земли нормальная компонента плотности тока из атмосферы должна обращаться в нуль, что в силу изотропии проводимости воздуха эквивалентно нулевой нормальной компоненте электрического поля, т, е, условию
д\/_
дп
= 0. (9)
т=Ее +И(в,ф)
В общем случае, когда проводимость под землей сравнима с проводимостью воздуха, необходимо решать единую задачу электропроводности, которая может быть записана как пара задач — для атмосферы и для самой Земли, с условиями непрерывности потенциала и нормальной компоненты плотности тока на поверхности раздела,
В настоящей работе рассматриваются краевая задача (5)-(7), хотя созданный алгоритм и его программная реализация позволяют решать задачи в слое, не являющемся сферическим, а также некоторые смешанные краевые задачи с условиями (8), (9) на разных участках границы.
3. Энергетический метод
Для задачи Дирихле (5)-(7) при естественных ограничениях на величину проводимости
0 < а1 < а < а2 < то, (10)
где а1, а2 — константы, справедлив принцип минимума функционала энергии [7], в соответствии с которым обобщенное решение краевой задачи может быть получено в результате минимизации функционала энергии
Ш(V) = У а^гав V)2 ¿П (11)
на множестве функций V, удовлетворяющих условиям (6), (7), Интеграл равен джоу-левой диссипации, т, е, выделению тепловой энергии, сопровождающему прохождение электрического тока по проводнику. Этим объясняется название функционала энергетическим, Линейная часть в функционале отсутствует, поскольку решается однородная задача. Квадратичная часть Ш(V) соответствует энергетическому скалярному произведению
а grad V grad и ¿П. (12)
Эта симметричная билинейная форма на рассматриваемом множестве функций положительно определена. Минимум функционала энергии существует, единствен и дает обобщенное решение задачи (5)-(7), которое при дополнительном предположении гладкости является классическим.
Энергетическая норма эквивалентна норме пространства Ш2(1)(П), поэтому обобщенное решение можно аппроксимировать кусочно-линейными функциями.
4. Вариационно-разностная схема
Для использования кусочно-линейных функций необходимо предварительно разбить область на тетраэдры. Считаем, что такое разбиение выполнено с соблюдением стандартных ограничений для нерегулярной сетки, которые обеспечивают эквивалентность Ь2 нормы кусочно-линейной функции и евклидовой нормы набора ее узловых значений Ут\
При аппроксимации энергетического пространства кусочно-линейными функциями получается вариационно-разностная схема. Уравнения для значений искомой функции V в узлах сетки получаются как условия минимума функционала энергии и представляют собой систему линейных алгебраических уравнений
Каждое из этих линейных алгебраических уравнений связывает значения V в данном узле и во всех соседних узлах. Матрица Ь получается симметричной, так как ее коэффициенты равны вторым производным функционала энергии, вычисленного для кусочно-линейной функции по узловым значениям этой функции, В силу первого неравенства (13) матрица является положительно определенной, поскольку кусочно-линейные функции принадлежат энергетическому пространству, в котором положительная определенность оператора доказана. Оценки собственных чисел этой матрицы только множителями а\ ъ о2 отличаются от оценок в случае уравнения Пуассона,
Вариационная формулировка задачи также позволяет обычным образом построить многосеточный метод. Использование для преобразования правых частей при укрупнении сетки оператора, сопряженного к оператору интерполяции, автоматически сохраняет симметрию и положительную определенность матриц на крупных сетках. Крупные сетки строятся путем выбрасывания каждого второго узла сетки по каждому направлению до тех пор, пока не останется всего один шаг сетки в одном из направлений, т, е, внутренних узлов не будет. На основной и вспомогательных сетках сглаженные решения строятся с помощью 5-10 итераций метода верхней релаксации с параметром т около 1.5,
5. Построение сетки
Удобно иметь сетку с такой же регулярной нумерацией узлов, как в параллелепипеде. Данные области будем называть криволинейными шестигранниками.
При моделировании атмосферы в целом приходится решать задачу в области, близкой к шаровому слою, поэтому сетка строится из лучей, проведенных из центра Земли через некоторые точки на сфере. Пример разбиения по высоте на 16 слоев показан на рис, 1, Сетка на сфере строится, например, проецированием равномерных прямоугольных сеток на гранях куба, имеющего тот же центр, что и Земля, В этом случае получается шесть криволинейных четырехугольников на поверхности Земли и соответственно шесть криволинейных шестигранников в атмосфере. Если Землю считать шаром, то все фигуры одинаковы. На рис, 2 показана часть такой сетки для полусферы, Приведена небольшая часть линий реальной сетки. Жирными линиями выделены
(13)
Чт * т
(14)
екая проекция
криволинейные четырехугольники. При представлении результатов используется стереографическая проекция полусферы на плоскость.
Построение структурированной сетки с предварительным разделением области на криволинейные четырехугольники или шестигранниками в областях, аналогичных кругу и сфере или шаровому слою и шару, позволяет избежать излишнего измельчения сетки и особенностей, возникающих около полюсов при использовании прямоугольных сеток в полярных или сферических координатах, а также сложностей нумерации соседних узлов нерегулярных сеток.
Дня использования кусочно-линейных функций каждую кубическую, в смысле нумерации ее узлов, ячейку сетки следует разбить на тетраэдры, что делается с предварительным определением некоторого центра ячейки и центров ее шести граней. Тогда каждая грань естественным образом делится на четыре треугольника и ячейка — на 24 тетраэдра. Значения кусочно-линейной функции в этих центрах получаются некоторой интерполяцией узловых значений, обычно как их среднее арифметическое.
6. Сходимость к точным решениям
При доказательство сходимости решения дискретной задачи к точному ограничимся областями П, являющимися многогранниками. Полагаем, что разбиение на тетраэдры выполнено с соблюдением стандартных условий дня кусочно-линейных восполнений на нерегулярной сетке, и обозначим через к максимальную длину отрезка в получившейся сетке.
Предположим, что существует решение задачи (5)-(7) — функция а, принадлежа-(2)
щая пространству Ш (П), Построим ее кусочно-линейное восполнение а., положив а.. = а в узлах сетки. Для погрешности аппроксимации справедлива оценка
II аи - а Ц(1)(П)< ск II а 11<)(П)> (15)
где константа с не зависит от шага сетки к [8].
Функционал энергии (11) может быть представлен в виде
ш(V) = [V - а, V - а] - [а, а],
(16)
где а — решение краевой задачи [7]. Возьмем в качестве V кусочно-линейное восполнение аъ функции а. При минимизации функционала энергии по всем возможным получим дискретное решение V?- Значение функционала энергии Ш (V?) ПРИ минимизации не возрастает, поэтому
[V0 - а, V? - а] < а - а, ан - а]. (17)
В силу эквивалентности энергетической нормы норме Ш2(1)(П) правая часть (17) оценивается сверху левой частью (15) с некоторым конечным множителем с1; зависящим от формы области и констант, ограничивающих проводимость в (10), но не зависящим от конкретной функции, нормы которой вычисляются. Получаем неравенство
([V? - а, V? - а])1/2 < С1сй II а ||^(2)(п), (18)
означающее пропорциональную к сходимость дискретных решений к точному, В силу эквивалентности энергетической нормы норме Ш2(1)(П) из этого неравенства следует такая же оценка сходимости в пространстве Ш2(1)(П),
Для получения фактической сходимости решим серию задач, имеющих точное решение.
7. Точные решения для сферически симметричной атмосферы
Когда область — шаровой слой и проводимость а зависит только от радиуса, уравнение (5) принимает вид
1 д ( 2 , \ а(г) д ( BV\ а(г) d2V . ч
--7Г rV гЬг " "Т^Т^ sin9^ " 2 2л д 2 = (19)
r2 д^ \ дг / r2 sin в дв \ дв j r2 sin2 в д^2
и его решение несложно получить методом разделения переменных. Представим V(r,e,^) в виде произведения F (r) на Y (в, и разделим уравнение на F (r)Y (в, ^)a(r)/r2, В результате получим
1 д ( ^ ^F(r) ' Г (7 (г)
F(r)a(r) dr \ dr
1 pe)
Y(в,^) Vsin вдв V дв y sin2 в %
Поскольку первое слагаемое зависит только от r, а второе — только от в, это равенство может быть выполнено лишь в случае, когда оба слагаемых — константы. Для Y (в, получается такое же уравнение, как и в случае уравнения Пуассона, и его частными решениями являются сферические функции, равные присоединенным полиномам Лежандра P™(cos в), умноженным на cos (m^) или на sin(m^), —n < m < n, Для сферических функций второе слагаемое (20) равно n(n + 1) [9]. Тогда для F(r) из (20)
получим обыкновенное дифференциальное уравнение
^ = + (21)
mn
Граничное условие (7) эквивалентно условию
Fn(RE) = 0, (22)
а условие (9) означает задание
Fn(Rup) = F°. (23)
Численное решение задачи (21)—(23) несложно. Решаем задачу Коши, стартуя от Fn(RE) = 0 с произвольным значением производной. Получаем некоторую функцию Fn(r) и, в частности, значение Fn(Rup), Не нарушая (21), (22), удовлетворяем (23) с помощью нормировки
Fn(r) = Fn(r)Fn0/Fn(Rup). (24)
При n = 0 получим сферически-симметричное решение, для которого функцию R(r) проще получить как интеграл
ад=<=/(?£(?)• <25>
где c — константа, которую используем, чтобы удовлетворить (23).
На рис. 3 показаны высотные зависимости Fn(r) при n = 0, 50, 100, 200. Интервал между ближайшими нулями полиномов Лежандра убывает примерно как 1/n, а им определяется горизонтальный масштаб решения. Случай n > 100 не рассматривается, поскольку горизонтальный масштаб становится менее 100 км, что не имеет смысла в силу сделанных в раздело 3 упрощений физической модели. Тем не менее мелкомасштабные решения для n = 500, 1000, 2000 на рисунке показаны, чтобы продемонстрировать их быстрое убывание с уменьшением высоты.
Построенные функции Fn(r) позволяют получить решение задачи (19), (22), (23) в виде ряда
V (r,0,p) = aoFo(r) +
<х / n \
+ Fn(r) anPn(cos 0) + Pnm(cos 0)(aw cos (m^) + bnm sin (m^)) , (26)
n=1 \ m=1
Рис. 3. Высотные зависимости частных решений ^П(г) при п = 0, 50,100,200, соответствующих полиномам Лежандра, и при п = 500,1000, 2000, соответствующих мелкомасштабным решениям
где коэффициенты anm, bnm получаются при разложении заданной функции V0(0, по сферическим гармоникам, В силу ортогональности такого базиса [9]
(2n + 1)(n — m)' Г
¿W = 1 0 , Д,-г^ / V0(e, V)Pu (cosd) cos (mp) dcos0 dip. (27)
2n(n + m)! J
Интегрирование производится по всей сфере 0 < 0 < п, 0 2п,ав выражении
для bnm только cos (m^) заменяется на sin (ш<^). Выражение для an получаем, полагая m = 0 и р0 = Pn-
Для тестирования алгоритма и программы ограничимся только частными решениями, не зависящими от т. е. полиномами Лежандра Fn(r)Pn(cos 0), поскольку иные m=0
8. Фактическая сходимость к точным решениям
Для тестирования программ и анализа точности получаемых решений были проведены тестовые расчеты для построенных в предыдущем разделе решений вида
V (r,0) = Fn(r)Pn(cos 0). (28)
Соответствующие значения потенциала изменяются с высотой от нуля до одного вольта. Использованы характерные значения n = 1,10,100,
Чтобы отличия от точных решений были видны, применена нарочито грубая сетка, содержащая всего восемь шагов по высоте. Это разбиение по высоте получается, если объединить пары слоев, показанных на рис, 1, Шаги по горизонтали подобраны достаточно мелкими, чтобы даже для самого мелкомасштабного решения, n = 100, не вносилась существенная дополнительная погрешность. Оказалось, что для такого разбиения целесообразно использовать шаг по горизонтали около 15 км,
Погрешность потенциала при n = 100 максимальна па высоте около 20 км, На рис, 4 показано распределение этой погрешности по горизонтали. Представленный квадрат 600 х 600 км является стереографической проекцией горизонтального сечения всей
Рис. 4. Распределение погрешности потенциала Уъ - V на горизонтальной поверхности на высоте 20 км; интервал между линиями уровня 0.005; штрих — отрицательные значения; сетка 40 х 40 х 8 п = 100; выделен квадрат 400 х 400 км
расчетной области. При этом на вертикальных границах значения потенциала задавались интерполяцией но вертикальным прямым значений с верхней границы. В качестве интерполяционной функции использовалось решение одномерной задачи, соответствующее (28) при п = 0. Если задать на этих границах точное решение, приграничные области с большой погрешностью исчезнут. При отходе от границы приграничная погрешность уменьшается почти экспоненциально, примерно в 300 раз на 100 км. Такой квадрат выделен па рис. 4. Малость ширины этих областей характеризует возможность проводить расчеты в небольших областях вместо расчетов сразу дня всей атмосферы. В частности, вместо показанной на рис. 2 совокупности криволинейных четырехугольников можно вести расчеты в каждом из них отдельно, предварительно расширив его всего на 200 км в двух горизонтальных направлениях.
Для рис. 4 использовалась грубая сетка, имеющая восемь слоев по высоте и 40 х 40 шагов но горизонтали. Отличие значений погрешности па соседних линиях уровня составляет 0.005, т. е. около половины процента от максимального значения потенциала. Максимальная погрешность наблюдается в центре области, где она достигает 3%.
При мельчепии сетки вдвое приграничная погрешность не меняется, а погрешность в основной части области убывает в 3.7 раза и при еще одном удвоении сетки — в 3.8 раза,
4
Дня решений с большими горизонтальными масштабами убывают обе составляющие погрешности, приграничная — поскольку вертикальная зависимость ближе к одномерному решению, в основной части области — поскольку получается больше точек па характерном масштабе. Так, при переходе от п = 100 к п = 50 первая по сравнению с
3.5
На рис. 5 приведены высотные распределения Ег (к) приближенных решений для (28) с параметром п = 100, полученные на разных сетках. Само точное решение не показано, так как в масштабе рисунка оно сливается с ломаной, соответствующей решению,
32
при мельчепии сетки. Штрих на рисунке — одномерное решение, которое получается
к, км
60 —-Х----
40 ---
\ 32\Л\8 20 - --
0 5 • 10"6 -Ег, В/м 10"5
Рис. 5. Высотные зависимости — Ег иа вертик али 9 = 0 да я п = 100 при максимальной разности потенциалов между землей и ионосферой, равной 1 В; три ломаных для сеток 40 х 40 х 8, 80 х 80 х 16 160 х 160 х 32 помечены количеством шагов по вертикали; штрих — одномерное решение
ч
1
\ 32^ \ ч
^¡х
из (28) при n = 0, В нижней части атмосферы оно является хорошим приближением для достаточно мелкомасштабного решения n = 100. При уменьшении n, т. е. при возрастании характерного горизонтального масштаба решения, соответствие еще более улучшается. Верхняя граница области адекватности одномерной аппроксимации поднимается от 20 км при n = 100 до 50 км при n = 10 и до верхней границы атмосферы n=0
На рис. 6 представлена скорость сходимости мпогосеточпых итераций дня задачи n = 100 при использовании сеток 40 х 40 х 8, 80 х 80 х 16 160 х 160 х 32. Узлы самой мелкой сетки по высоте распределены по тому же правилу, что и па рис. 1 дня сетки 80 х 80 х 16
же область. В качестве нормы невязки па левом фрагменте рисунка используется сумма модулей невязок алгебраических уравнений во всех узлах сетки, характеризующая дисбаланс токов и выражающаяся в амперах. Отметим, что полный ток, проходящий через расчетную область, составляет величину порядка 0.01 А,
Начальное приближение строилось интерполяцией с помощью одномерных решений. Линейная по высоте начальная интерполяция увеличивает начальную норму невязки па порядок. Если начинать итерации с нулевого потенциала внутри области, то начальная норма невязки будет примерно в сто раз больше. Иначе говоря, удачный выбор начального приближения позволяет сэкономить пару итераций. Дня самой манкой сетки па рис. 6, а топкой ломаной дополнительно показана сходимость при нулевом начальном приближении.
Как видим, дня самой грубой сетки скорость сходимости несколько выше, а дня более мелких, как и предсказывается теорией мпогосеточпых методов, — становится независимой от числа узлов сетки. Дня мпогосеточпых итераций вспомогательные крупные
0
log || Я -5
-10 -15
-20
Рис. 6. Сходимость многосеточных итераций для задачи п = 100; по горизонтали — номер итерации; три ломаных для сеток 40 х 40 х 8 80 х 80 х 16 160 х 160 х 32 помечены количеством шагов по вертикали: а норма невязки; для последней сетки тонкой ломаной показана
13
горизонтали сходимость итераций Зейделя при нулевом (кривая 0) и одномерном (кривая 1) начальном приближении; б — норма погрешности потенциала Уь — У'- тонкими ломаными показано отличие У^ после данной итерации от того У^, которое установится после 20 итераций
сетки строились путем выбрасывания каждого второго узела сетки по каждому направлению до тех пор, пока не исчезнут все внутренние узлы. Соответственно для основных сеток 40 х 40 х 8, 80 х 80 х 16 160 х 160 х 32 использовались три, четыре и пять
10
лакеации с параметром т, оптимальное значение которого меняется от т = 1.25 для самой крупной из перечисленных сеток до т = 1.45 и 1.5 для более мелких. Уменьшение вдвое количества этих итераций на каждой сетке ведет примерно к удвоению количества необходимых многосеточных итераций, что практически сохраняет время расчетов.
На рис, 6, а штрихом приведена сходимость итераций Зейделя для сетки 160 х 160 х 32 при нулевом и одномерном начальном приближении. Каждая рассматриваемая мно-
13
еравнения эффективности итераций штриховые ломаные, характеризующие итерации
13
сетке примерно в 6 раз эффективнее итераций Зейделя, С ростом числа узлов выигрыш увеличивается.
Как следует из рис, 6, о, за пять итераций норма невязки убывает до минимально возможного значения, которое обусловлено конечной точностью машинных чисел, В расчетах использовалась арифметика двойной точности, т, е, восьмибайтовые чис-16
Разноеть между приближенным Т4 и точным V решениями можно характеризовать ее нормой в пространстве ^(П). Разделив эту норму на объем области |П|, получим более наглядное среднее значение модуля функции Vh — V:
Приближенное значение этого интеграла вычисляется, полагая функции в каждой ячейке сетки постоянными и равными узловому значению. Тогда интеграл превращается в сумму разностей узловых значений с весами, равными объемам ячеек. Чтобы исключить приграничную погрешность, можно рассматривать — V только в некоторой подобласти, проведя расчеты в более широкой области. Например, при рассматриваемых здесь расчетах для п = 100 в области 600 х 600 км в подобласти 400 х 400 км получаем значения нормы (29) 0.00759, 0.00200 и 0.00052, Уменьшение погрешности примерно вчетверо при мельчении сетки вдвое демонстрирует сходимость приближенных решений к точному, близкую к квадратичной. Отметим, что сама функция V имеет значение порядка единицы.
На рис, 6, б"жирными ломаными показано убывание нормы погрешности потенциала Vh — V в процессе итераций, тонкими ломаными — отличие Vh после данной итерации от того Уь которое установится после 20 итераций. Видно, что сходимость процесса решения системы линейных алгебраических уравнений (14) в этой норме продолжается еще несколько итераций после того, как норма невязки уже достигла минимального значения. После ряда первых итераций погрешность приближенного решения Vh — Vyжe не уменьшается, поскольку она обусловлена погрешностью аппроксимации.
Рисунок 6, б демонстрирует, что для самой грубой сетки итерации практически не уменьшают норму Vh — V, которая остается той же, что была сразу после построения начального приближения. Это объясняется тем, что разность между точным и приближенным решениями V — ^ < 0.05 и разность между одномерным и точным
(29)
решениями V — У1-в < 0.08 различаются всего в полтора раза, что мало заметно в масштабе рис, 6, б. Однако малое изменение нормы погрешности потенциала в процессе итераций не означает, что итерации не нужны. Достаточно обратить внимание па вертикальную компоненту плотности тока ]г. Начальное приближение, которое строится как одномерное решение, соответствует нулевому полиному Лежапдра и имеет практически постоянную по высоте плотность вертикального тока jr = 10-17 А/м2. В верхней атмосфере отличие такой jr от рассматриваемого точного решения, соответствующего сотому полиному Лежапдра, достигает четырех порядков и почти устраняется итерациями.
Дня достижения практически той же точности всех представленных параметров, что и после большого количества итераций, достаточно двух мпогосеточпых итераций. На обычном персональном компьютере с частотой 2 ГГц они занимают время около 15 с для сетки 160 х 160 х 32, имеющей около миллиона узлов. Все массивы еще входят в оперативную память, которая используемым Фортраном ограничена 270 МБ. Гораздо больше времени, около 100с, занимает вычисление коэффициентов матрицы (14).
9. Расчет атмосферных электрических полей для спокойных геомагнитных условий
Распределения потенциала в ионосфере дня различных геомагнитных условий приводятся, например, в работе |10|, В представленных расчетах используется модельное распределение дня спокойных геомагнитных условий, полученное в модели |11|, Линии
5
Расчеты выполнены дня Северного полушария па сетке, показанной па рис. 2, по содержащей 513 х 513 х 17 узлов в центральном криволинейном шестиграннике и вдвое меньше — в остальных четырех. Так как основные ноля сосредоточены в высоких и средних широтах, то получаются те же результаты, если дня покрытия интересующей
Рис. 7. Распределение Ег на высоте 22 км жирные линии; разность значении Ег между соседними линиями уровня составляет 13.5 мВ/м; одномерное решение показано тонкими линиями, которые совпадают с эквипотенциалями в ионосфере с интервалом 5 кВ
нас области ограничиться одним центральным немного расширенным шестигранником. Поскольку шаги сетки соответствуют выбранным в тестовых расчетах предыдущего раздела, сходимость итераций и сходимость по шагу сетки не отличаются от тестовых. Сетка 513 х 513 х 17 содержит около 4.5 млн узлов. При независимых расчетах в каждой
100
етея в оперативной памяти 270 МБ, используемой имеющимся у авторов транслятором Фортрана, Тогда реальное время работы компьютера примерно равно суммарному вре-
3
процессора, как и размеры массивов, увеличивается вчетверо, а реальное время рабо-
1000
12 500 70
ней памяти имеет место шестикратное замедление, С увеличением числа узлов сетки аналогичное замедление происходит и на компьютере с большей памятью. Поэтому целесообразно вести независимые расчеты в подобластях. Результаты расчетов, как и в описанных выше тестах, не меняются.
На рис, 7 жирными линиями показано распределение полученной вертикальной компоненты электрического поля Ег на высоте 22 км. Интервалы между линиями уровня Ег и потенциала подобраны так, чтобы для одномерного решения эти линии совпадали, В нижней атмосфере это приближение хорошо работает, поскольку до высоты около 50 км разрезы Ег остаются похожими представленным на рис, 7, Одномерное решение соответствует постоянству плотности вертикального тока, точнее, ее изменению с высотой как 1/г2, которым можно пренебречь. Поэтому на уровне земли Ег получается умножением на величину а(20 км)/а(0). Максимальное значение, около 30В/м, является заметным изменением наблюдаемого среднего поля 130 В/м,
При более высокой геомагнитной активности разности потенциала в ионосфере мо-100
60
чины наблюдаемого поля.
Существенным свойством исходного распределения потенциала на верхней границе атмосферы является его гладкость. Если использовать исходное распределение потенциала на достаточно грубой сетке и интерполировать его на расчетную сетку, то результаты для верхней атмосферы в существенной степени определяются выбранной интерполяцией. При линейной интерполяции имевшиеся скачки производных на границах ячеек исходной сетки проявляются и у потенциала на расчетной мелкой сетке. Картина эквипотенциалей выглядит довольно гладкой (см, тонкие линии на рис, 7), но в разложении такой функции есть высокие гармоники. Как видно из рис, 3, увеличение порядка полинома Лежандра дает быстрый рост Ег, т. е, вертикальной производной потенциала. Решение существенно меняется. Например, в горизонтальных сечениях вертикальной компоненты электрического поля Ег в верхней половине атмосферы отчетливо проявляются все ребра исходной грубой сетки, которая использовались для расчетов ионосферного поля в [11]. Аналогичный приведенному на рис, 7 вид сохраняется до тем больших высот, чем более гладким берется исходное распределение потенциала на верхней границе атмосферы. Если провести сто итераций сглаживания исходного распределения потенциала, каждый раз заменяя его значение в узлах на средние
60
показанное на рис, 8, При таком сглаживании увеличивается радиус кривизны эквипотенциалей в местах их наибольшего изгиба. Участки линий при этом смещаются на
Рис. 8. Распределение Ег на высоте 60 км (жирные линии); разность значении Ег между соседними линиями уровня составляет 69мкВ/м; одномерное решение показано тонкими линиями,
5
расстояния порядка расстояний между линиями. Подчеркнем, что рис. 7, 8 относятся к разным решениям, отличающимся сглаживанием исходных данных.
Если рассматривать линии уровня Ег на фиксированной большой высоте, где картина линий уровня совершенно неудовлетворительна, то со сглаживанием исходного потенциала исчезают мелкомасштабные погрешности и картина становится гладкой, по сильно отличается от одномерной. Как видно из рис. 8, ноле имеет иное пространственное распределение и его напряженность в областях экстремумов удвоена по сравнению с одномерным решением, показанным топкими линиями. Следует отметить, что грубая расчетная сетка тоже увеличивает долю верхней атмосферы, в которой сказываются погрешности задания потенциала па верхней границе.
Таким образом, перед расчетами глобального ноля необходимо выделить крупномасштабную, или, иначе говоря, гладкую, часть. Отброшенную мелкомасштабную часть исходного распределения потенциала па верхней границе атмосферы, если она пе является просто погрешностью представления функции, целесообразно рассмотреть отдельно, проведя расчеты в небольших подобластях.
10. Параллельные вычисления
Представленный в разделе 5 способ построения сетки позволяет дополнительно разрезать любой криволинейный шестигранник, В описанном в разделе 4 алгоритме основной объем вычислений требуется при проведении итераций метода верхней релаксации па основной и вспомогательных сетках. Их можно проводить независимо внутри каждого шестигранника с обменом сравнительно небольшой информацией около границ, в силу чего предложенный алгоритм может быть эффективно распараллелен. Это позволяет практически линейно сокращать время расчетов в зависимости от числа процессоров, а также использовать более мелкую сетку дня повышения точности решения.
Проверка эффективности параллельных вычислений проводилась для области, описанной в разделе 5, дополнительно разделенной по горизонтали. Центральный шестигранник был разделен на четыре части, а каждый боковой — на две, В каждом из получившихся 12 шестигранников построена сетка 256 х 256 х 16, Для решения всей задачи на одном компьютере потребовалось бы около трех гигабайт оперативной памяти, В нашем комплексе программ требуется большой объем памяти, поскольку он рассчитан на достаточно произвольные сетки, когда уравнения системы линейных алгебраических уравнений вариационно-разностной схемы связывают значения потенциала в узле
26
14
10
48
процессоров с раздельной памятью по одному гигабайту на процессор.
Количество процессоров 2 3 4 6 12 Время расчетов, с 545 155 116 80 41
В случае использования двух процессоров видно значительное увеличение времени расчетов, связанное с физическими ограничениями применяемого компьютера. Не полностью вместившаяся в оперативную память компьютера область частично сохранялась в своп-файле. При использовании трех и более процессоров эта проблема не возникала.
Как видно из приведенных данных, время расчетов убывает практически обратно пропорционально числу используемых процессоров. Необходимо отметить, что для равномерной загрузки процессоров количество примерно равных шестигранников должно быть кратно количеству процессоров.
Выводы
Построен эффективный вычислительный алгоритм для решения глобальной задачи электропроводности атмосферы.
Предложено проводить не единый расчет, а серию расчетов в небольших подобластях, что существенно уменьшает время вычислений за счет исключения обращений к внешней памяти компьютера. Показано, что для моделирования локальных полей и то-
100
ях и, чтобы не возмущать решение в интересующей нас области атмосферы, поставить на вертикальных границах условия Дирихле с помощью одномерных решений. Преимущество по сравнению с методом декомпозиции области здесь состоит в отсутствии дополнительных итераций, требующих неоднократных решений задач в подобластях.
Показано, что электрическое поле вблизи поверхности Земли можно строить как решение одномерной по высоте задачи, но чтобы получить поля и токи в верхней атмосфере, необходимо решать трехмерную задачу.
Представленный алгоритм легко распараллеливается и может выполняться на многопроцессорных компьютерах, что позволяет в силу уменьшения шага сетки сокращать время расчетов и получать более точное решение.
Список литературы
fl] Parrot М. Use of satellites to detect seismo-electromagnetic effects // Adv. Space Res. 1995. Vol. 15. P. (11)27—(11)35.
[2] Molchanov O., Fedorov E., Schekotov A. et al. Lithosphere-atmosphere-ionosphere coupling as governing mechanism for preseismic short-term events in atmosphere and ionosphere // Nat. Hazards Earth Svs. Sci. 2004. Vol. 4. P. 757-767.
[31 Denisenko V.V., Boudjada M.Y., Horn M. et al. Ionospheric conductivity effects on electrostatic field penetration into the ionosphere // Ibid. 2008. Vol. 8. P. 1009-1017.
[4] Molchanov O., Hayakawa M. Seismo-Electromagnetics and Related Phenomena: History and Latest Results. Tokyo: TERRAPUB, 2008. 190 p.
[5] Denisenko V.V., Biernat H.K., Mezentsev A.V. et al. Modification of conductivity due to acceleration of the ionospheric medium // Ann. Geophvs. 2008. Vol. 26. P. 2111-2130.
[6] Денисенко В.В. Энергетический метод для трехмерных эллиптических уравнений с несимметричными тензорными коэффициентами // Сибирский мат. журн. 1997. Т. 38, № 6. С. 1267-1281.
[7] Михлин С.Г. Вариационные методы в математической физике. М.: Гостехиздат, 1957. 378 с.
[8] Оганесян Л.А., руховец Л.А. Вариационно-разностные методы решения эллиптических уравнений. Ереван: Изд-во АН Арм. ССР, 1979. 235 с.
[9] Янке Е., Эмде Ф., Леш Ф. Специальные функции. Формулы, графики, таблицы. М.: Физматгиз, 1968. 344 с.
[10] Weimer D.R. Improved ionospheric electrodvnamic models and applications to calculating Joule heating rates //J. Geophvs. Res. 2005. Vol. 110. P. A05306.
[11] Denisenko V.V., Zamay S.S. Electric field in the equatorial ionosphere // Planetary Space Sci. 1992. Vol. 40, No. 7. P. 941-952.
Поступила в редакцию 9 октября 2009 г.
