авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 || 3 | 4 |

«ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ УЧРЕЖДЕНИЕ НАУКИ ИНСТИТУТ ВЫЧИСЛИТЕЛЬНОГО МОДЕЛИРОВАНИЯ СИБИРСКОГО ОТДЕЛЕНИЯ РОССИЙСКОЙ АКАДЕМИИ НАУК ...»

-- [ Страница 2 ] --

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

Есть иной подход, базирующийся на аппроксимации производных конеч ными разностями [2]. В его рамках обычно обосновывается более сильная, поточечная сходимость, которая получается как следствие аппроксимации и устойчивости.

В рассматриваемом частном случае при подстановке рядов Тейлора ви да (2.59-2.61) в уравнение (2.57) все члены с нечетными степенями h и без множителя h взаимно сокращаются. После традиционного деления на h3 по лучаем 2V 2V 2V qnml + h2 O(1) = 3, (2.63) + + x2 y 2 z 2 h где остаточный член содержит некоторую комбинацию четвертых производ ных функции V (x, y, z), вычисленных в некоторых точках куба |x| h, |y| h, |z| h. В предположении ограниченности всех четвертых произ водных эта величина порядка 1, как это записано в (2.63).

В разностных схемах обычно правая часть вычисляется просто как правая часть дифференциального уравнения в данном узле сетки, то есть (2.64) qnml = q(x, y, z), h и тогда погрешность ее аппроксимации нулевая.

В построенной нами вариационно-разностной схеме (2.38) qnml получает ся перераспределением qi из вспомогательных узлов по формулам (2.37). В рассматриваемом частном случае, когда ячейки сетки являются кубами, цен тры ячеек (точка 15 на рис. 2.5) и граней (точка 9 на рис. 2.5) определяются соображениями симметрии как центры куба и квадрата, коэффициенты ин терполяции aik в (2.31) равны 1/8 и 1/4. Поэтому в силу (2.37) qi из центра ячейки поровну перераспределяется в qnml восьми вершин ячейки, а из цен тра грани – поровну в четыре вершины данного квадрата. Для того, чтобы вычислить сами qi, в силу (2.18) надо провести интегрирование (2.28).

Для численной реализации этого интегрирования есть много подходов.

Например, предварительно заменив функцию q(x, y, z) функцией, постоян ной в каждом тетраэдре сетки, получаем последнее выражение (2.28). Вы числение правых частей нашей вариационно-разностной схемы является гро моздкой процедурой, если буквально следовать теории.

На практике мы в некоторых случаях используем простейшее приближение (2.64). Поскольку в задачах электропроводности q(x, y, z) есть минус дивергенция сторонних токов, естественно приближенно вычислять ток между соседними ячейками сетки, например, через полуцелую грань n + 1/2, m, l, и вносить этот ток с плюсом в qn1,ml и с минусом в qnml. При этом точно выполняется закон сохранения заряда, проинтегрированный по всей области. В случае задачи Неймана это свойство является необходимым условием разрешимости зада чи. Отметим, что в описанных в диссертации моделях сторонние токи не рассматриваются, и поэтому q(x, y, z) 0, хотя в следующем разделе мы опишем процедуру вычисления qnml при наличии сторонних токов.

Значения Vnml в граничных узлах мы задаем точно нулевыми и поэто му погрешность аппроксимации в смысле теории разностных схем [2] равна нулю.

Таким образом, в частном случае, когда проводимость постоянна, q(x, y, z) = 0, и сетка составлена из одинаковых кубических ячеек, постро енная вариационно-разностная схема (2.38) принимает вид (2.57) с нулевой правой частью и является разностной схемой второго порядка аппроксимации для задачи электропроводности (1.17, 2.1, 2.3).

Еще одно важное свойство получившейся разностной схемы – устойчи вость, то есть оценка сверху max|Vnml | через max|qnml /h3 |, установлена нера венством (2.46). Из доказанных аппроксимации и устойчивости следует схо димость сеточных значений к точному решению в узлах сетки [2].

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

2.9. Дискретная модель.

Вариационно-разностной схеме (2.38) можно придать наглядную форму за кона сохранения заряда в интергральной форме. Начнем с выражения функ ционала энергии (2.18). Каждое из уравнений (2.29) получается дифферен цированием W (V h ) по параметру Vi :

I (2.65) Vj grad vi · grad vj d = (viq + vi · )d.

j= Векторы grad vi отличны от нуля только в тетраэдрах, прилегающих к рассматриваему i-тому узлу. Рассмотрим один такой тетраэдр.

С учетом (2.22), внося суммирование под знак интеграла, преобразуем левую часть:

I Hi Hi Hi h h (2.66) 2 · grad ( Vj vj )d = 2 · grad V d = · j d, Hi Hi Hi j= где через jh обозначена соответствующая потенциалу V h плотность тока.

При проводимости, постоянной внутри каждого тетраэдра, вектор jh тоже постоянен, и значит интеграл (2.66) равен tetr S Hi Hi · jh · jh = jh · S, (2.67) = Hi2 Hi 3 где S - вектор площади треугольника, на который опущена высота Hi.

i Рис. 2.6. Разделение граней тетраэдра на треугольники.

Последнее скалярное произведение есть сила тока через треугольник S. В силу постоянства jh этот ток имеет нулевую дивергенцию, и значит jh · S ра вен со знаком минус суммарному току через остальные три грани тетраэдра.

Разделим каждую из них на три равные по площади части, взяв середины отрезков и некоторые центры треугольников, как это показано на рис. 2.6.

Теперь последнее выражение (2.67) можно трактовать как взятый с обрат ным знаком суммарный ток через три получившихся малых треугольника, прилегающих к узлу i. А он, в силу бездивергентности jh, равен току через произвольную поверхность, проведенную внутри тетраэдра через показанную пунктирами ломанную, окружающую i-й узел. Удобно взять некоторый центр тетраэдра и соединить его плоскими треугольниками со всеми отрезками, по казанными пунктиром на рис. 2.6. Центр выбираем так, что равны объемы всех четырех частей, на которые оказывается разделен тетраэдр. Теперь по следнее выражение (2.67) есть ток через построенную поверхность, которая ограничивает прилегающую к i-тому узлу четверть тетраэдра. Суммируя по всем прилегающим к i-тому узлу тетраэдрам, получаем ток из выделенной окрестности этого узла. Причем вся область оказалась разбита на такие окрестности узлов.

Аналогичные рассуждения для последнего интеграла (2.65), содержащего вспомогательную функцию, позволяют трактовать этот интеграл как соот ветствующую потенциалу часть тока, вытекающую из той же окрестности i-того узла. То же можно проделать и для оставшейся заданной функции q, вспомнив, что она есть дивергенция стороннего тока (1.4). Поэтому справед лива цепочка равенств:

vidiv jext d = div(vi jext )d + (2.68) qvi d = jext · grad vi d.

Последний интерграл имеет тот же вид, который позволил получить вы ражение (2.67) и придать ему смысл суммарного тока из окрестности i-того узла, только теперь речь идет о стороннем токе. Предпоследний интеграл равен интегралу по границе (2.69) vi jext dS, и равен нулю, поскольку все функции vi обращаются в нуль в узлах границы и на всех граничных треугольниках сетки.

Таким образом, мы придали уравнениям вариационно-разностной схемы (2.38) смысл закона сохранения заряда, проинтегрированного по окрестно сти каждого узла сетки. Естественно, речь идет о полном токе, равном сум ме стороннего и тока проводимости, который мы из чисто математических соображений разделили на две части, соответствующих представлению элек трического потенциала V в виде суммы двух функций V = V + (2.12).

Переход от алгебраических уравнений (2.29) к (2.38) происходит при ис ключении вспомогательных узлов с помощью интерполяции (2.31). С этими же коэффициентами aik соответствующие уравнения (2.29) перераспределя ются в остающиеся. Поэтому получаются законы сохранения заряда для бо лее сложно устроенных окрестностей узлов сетки. Например, если сетка со ставлена из кубических ячеек, то к построенной выше окрестности i-того узла добавляются четверти окрестностей узлов типа 9 на рис. 2.5 и 1/8 окрестно стей узлов типа 15 на том же рисунке.

2.10. Сходимость к точным решениям.

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

Предположим, что существует решение задачи (1.17, 2.1, 2.2) – функция (2) a, принадлежащая пространству W2 (). Построим ее кусочно-линейное вос полнение ah, положив ah = a в узлах сетки. Для погрешности аппроксимации справедлива оценка ah a W (1) () c3 h a W (2) (), (2.70) 2 где константа c3 не зависит от шага сетки h [32].

Функционал энергии (2.13) может быть представлен в виде (2.71) W (V ) = [V a, V a] [a, a], где a – решение краевой задачи [28]. Возьмем в качестве V кусочно-линейное восполнение ah функции a. При минимизации функционала энергии по всем возможным V h получаем дискретное решение V0h. Значение функционала энергии W (V0h ) при минимизации не возрастает, поэтому [V0h a, V0h a] [ah a, ah a]. (2.72) (1) В силу эквивалентности энергетической нормы норме W2 () правая часть (2.72) оценивается сверху левой частью (2.70) с некоторым конечным множителем c4, зависящим от формы области и констант, ограничивающих проводимость в (2.5), но не зависящим от конкретной функции, нормы кото рой вычисляются. Получаем неравенство ([V0h a, V0h a])1/2 c4 c3 h a W (2) (), (2.73) означающее пропорциональную h сходимость дискретных решений к точно (1) му. В силу эквивалентности энергетической нормы норме W2 () из этого (1) неравенства следует такая же оценка сходимости в пространстве W2 ().

Для выяснения фактической сходимости решим серию задач, имеющих точное решение.

2.11. Точные решения для сферически симметричной атмосферы.

Когда проводимость зависит только от радиуса, уравнение (1.17) принимает вид (r) 2V 1 V (r) V (2.74) 2 r (r) 2 sin 2 2 = 0.

r sin r r r r sin Когда область является шаровым слоем, решение краевой за дачи несложно получить методом разделения переменных. Представив V (r,, ) в виде произведения F (r) на Y (, ), и разделив уравнение на F (r)Y (, )(r)/r2, получаем 1 F (r) r2 (r) F (r)(r) r r 1 2Y (, ) 1 1 Y (, ) (2.75) sin + = 0.

sin2 Y (, ) sin Поскольку первое слагаемое зависит только от r, а второе – только от,, это равенство может быть выполнено, только когда оба слагаемых – константы. Для Y (, ) получается такое же уравнение, как и в случае урав нения Пуассона, и его частными решениями являются сферические функ m ции, равные присоединенным полиномам Лежандра Pn (cos ), умноженным на cos (m) или на sin (m), n m n. Для сферических функций вто рое слагаемое (2.75) равно n(n + 1) [47]. Тогда для F (r) из (2.75) получаем обыкновенное дифференциальное уравнение d dFn(r) r2 (r) (2.76) = n(n + 1)Fn(r)(r), dr dr которое не зависит от индекса m, и поэтому нумеруется только индексом n.

Граничное условие (2.2) эквивалентно условию (2.77) Fn (RE ) = 0, а условие (2.4) означает задание (2.78) Fn (Rup) = Fn.

Используемый нами метод численного решения задачи (2.76-2.78) подроб но описан в разделе 4.4. Сначала решаем задачу Коши, стартуя от Fn (RE ) = с единичным значением производной. Получаем некоторую функцию Fn (r), и в частности, значение Fn (Rup). Не нарушая (2.76, 2.77), удовлетворяем (2.78) с помощью нормировки (2.79) Fn (r) = Fn (r)Fn /Fn(Rup).

При n = 0 получаем сферически симметричное решение, для которого функцию F (r) проще получить как интеграл dr (2.80) F0 (r) = c, (r)2 (r) где c – константа, которую используем, чтобы удовлетворить (2.78).

На Рис. 2.7 показаны высотные зависимости Fn (r) при n = 0, 50, 100, 200.

Интервал между ближайшими нулями полиномов Лежандра убывает при мерно как 1/n, а им определяется горизонтальный масштаб решения. Мы не рассматриваем n 200, поскольку горизонтальный масштаб становит ся менее 100 км, что не имеет смысла в силу сделанных в разделе 3 упро щений физической модели. Тем не менее, мелкомасштабные решения, n = 500, 1000, 2000, тоже показаны на этом рисунке, чтобы продемонстрировать их быстрое убывание с уменьшением высоты.

Построенные функции Fn (r) позволяют получить решение задачи (2.74, 2.77, 2.78) в виде ряда V (r,, ) = a0 F0(r) + Fn (r) · n= n m (2.81) · an Pn (cos ) + Pn (cos )(anm cos (m) + bnm sin (m)), m= где коэффициенты anm, bnm получаются при разложении заданной функции V0 (, ) по сферическим гармоникам. В силу ортогональности такого бази са приведенное, например, в справочнике [47] выражение для интеграла от m квадрата Pn (6.6) позволяет получить (2n + 1)(n m)! m (2.82) anm = V0 (, )Pn (cos ) cos (m) d cos d.

2(n + m)!

h [км] 200 0.75 V 1. 0.00 0.25 0. Рис. 2.7. Высотные зависимости частных решений Fn (r), соответствующих полиномам Лежандра, при n = 0, 50, 100, 200, а также при n = 500, 1000, 2000, соответствующих мелкомасштабным решениям.

Интегрирование производится по всей сфере 0, 0 2, а в выражении для bnm только заменяется cos (m) на sin (m). Выражение для an получаем, полагая m = 0, и считая Pn = Pn.

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

Для тестирования алгоритма и программы ограничимся только част ными решениями, не зависящими от, то есть полиномами Лежандра Fn(r)Pn (cos ), поскольку иные решения, m = 0, получаются из них пово ротами.

2.12. Фактическая сходимость к точным решениям.

Для тестирования программ и анализа точности получаемых решений были проведены тестовые расчеты для построенных в предыдущем разделе реше ний (2.83) V (r, ) = Fn(r)Pn (cos ), изменяющихся с высотой от нуля до 1 вольта. Использованы характерные значения n = 1, 10, 100.

Использована нарочито грубая сетка, содержащая всего 8 шагов по высо те, чтобы отличия от точных решений были видны. Это разбиение по высоте получается, если объединить пары слоев, показанных на Рис. 2.2. Шаги по горизонтали подобраны достаточно мелкими, чтобы даже для самого мелко масштабного решения, n = 100, не вносилась существенная дополнительная погрешность. Оказалось, что для такого разбиения по высоте целесообразно использовать шаг по горизонтали около 15 км.

Рис. 2.8. Распределение погрешности потенциала Vh V на горизонтальной поверхности на высоте 20 км. Между линиями уровня 0.005, пунктир – отрицательные значения. Сетка 40 40 8, n = 100. Выделен квадрат 400 км 400 км.

Погрешность потенциала при n = 100 максимальна на высоте около км. На Рис. 2.8 показано распределение этой погрешности по горизонтали.

Показанный квадрат 600 км 600 км является стереографической проекцией горизонтального сечения всей расчетной области. При этом на вертикальных границах значения потенциала задавались интерполяцией по вертикальным прямым значений с верхней границы. В качестве интерполяционной функ ции использовалось решение одномерной задачи, соответствующее (2.83) при n = 0. Если задать на этих границах точное решение, приграничные обла сти с большой погрешностью исчезнут. При отходе от границы приграничная погрешность уменьшается почти экспоненциально, примерно в триста раз на 100 км. Такой квадрат выделен на Рис. 2.8. Малость ширины этих областей характеризует возможность проводить расчеты в небольших областях вместо расчетов для всей атмосферы сразу. В частности, вместо показанной на Рис.

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

Для Рис. 2.8 использовалась грубая сетка, имеющая 8 слоев по высоте и 40 40 шагов по горизонтали. Отличие значений погрешности на соседних линиях уровня равно 0.005, то есть около половины процента от максималь ного значения потенциала. Максимальна погрешность в центре области, где она достигает 3 процентов.

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

Для решений с большими горизонтальными масштабами убывают обе со ставляющие погрешности, приграничная – поскольку вертикальная зависи мость ближе к одномерному решению, в основной части области – поскольку получается больше точек на характерном масштабе. Так, при переходе от n = 100 к n = 50 первая уменьшается в 3.5 раза по сравнению с Рис. 2.8, вторая – втрое.

h [км] 32 5 · 106 Er, [В/м] Рис. 2.9. Высотные зависимости Er, [В/м] на вертикали = 0 для n = при максимальной разности потенциалов между землей и ионосферой, равной 1 [В]. Три ломаных для сеток 40 40 8, 80 80 16, 160 160 помечены количеством шагов по вертикали. Пунктир – одномерное решение.

На Рис. 2.9 приведены высотные распределения Er (r) приближенных ре шений для (2.83) с параметром n = 100, полученные на разных сетках. Само точное решение не показано, поскольку в масштабе рисунка оно сливается с ломаной, соответствующей решению, полученному на сетке с 32 слоями по высоте. Видна сходимость к точному решению при мельчении сетки.

Пунктиром на этом рисунке показано одномерное решение, которое полу чается из (2.83) при n = 0. Видно, что в нижней части атмосферы оно яв ляется хорошим приближением для достаточно мелкомасштабного решения n = 100. При уменьшении n, то есть при возрастании характерного горизон тального масштаба решения, соответствие еще улучшается. Верхняя граница области адекватности одномерной аппроксимации поднимается от 20 [км] при n = 100 до 50 [км] при n = 10, и, естественно, до верхней границы атмосферы при n = 0.

0 log H log Vh V 0 5 10 15 8 20 0 5 10 0 5 Рис. 2.10. Сходимость многосеточных итераций для задачи n = 100. По горизонтали – номер итерации. Три ломаных для сеток 40 40 8, 80 80 16, 160 160 32 помечены количеством шагов по вертикали. а – норма невязки. Для последней сетки тонкой ломаной дополнительно показана сходимость при нулевом начальном приближении, а также пунктиром в сжатом в 13 раз масштабом по горизонтали – сходимость итераций метода последовательной верхней релаксации при нулевом – кривая 0, и одномерном – кривая 1, начальных приближениях. б – норма погрешности потенциала Vh V. Тонкими ломаными показано отличие Vh после данной итерации от того Vh, которое установится после 20 итераций.

На Рис. 2.10 показана скорость сходимости многосеточных итераций для задачи n = 100 при использовании сеток 40 40 8, 80 80 16, 160 160 32.

Узлы самой мелкой из них по высоте распределены по тому же правилу, что показано на Рис. 2.2 для сетки 80 80 16. По горизонтали они все примерно равномерны и покрывают одну и ту же область.

В качестве нормы невязки на левом фрагменте используется сумма мо дулей невязок алгебраических уравнений во всех узлах сетки. Она характе ризует дисбаланс токов и выражается в амперах. Отметим, что полный ток, проходящий через расчетную область порядка 0.01 [А].

Начальное приближение строилось интерполяцией с помощью одномер ных решений. Линейная по высоте начальная интерполяция увеличивает на чальную норму невязки на порядок. Если начинать итерации с нулевого по тенциала внутри области, начальная норма невязки была бы примерно в сто раз больше. Как видим, удачный выбор начального приближения позволяет сэкономить пару итераций. Для самой мелкой сетки на Рис. 2.10 а тонкой ломаной дополнительно показана сходимость при нулевом начальном при ближении.

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

Соответственно, для основных сеток 40 40 8, 80 80 16, 160 160 использовались 3, 4, 5 вспомогательных сеток. На каждой сетке выполняется 10 итераций метода верхней релаксации с параметром, оптимальное зна чение которого меняется от = 1.25 для самой крупной из перечисленных сеток до 1.45 и 1.5 для более мелких. Уменьшение вдвое количества этих ите раций на каждой сетке ведет примерно к удвоению количества необходимых многосеточных итераций, что практически сохраняет время расчетов.

На Рис. 2.10 а пунктиром показана сходимость итераций метода после довательной верхней релаксации для сетки 160 160 32 при нулевом и од номерном начальных приближениях. Каждая наша многосеточная итерация по затратам эквивалентна 13 итерациям метода последовательной верхней релаксации. Отметим, что при значении релаксационного параметра = метод последовательной верхней релаксации становится методом Зейделя, который также называется методом Гаусса-Зейделя. Для удобства сравнения эффективности итераций пунктирные ломаные, характеризующие итерации метода последовательной верхней релаксации, сжаты в 13 раз по горизонта ли. Как видим, многосеточные итерации на этой сетке примерно вшестеро эффективнее. С ростом числа узлов выигрыш увеличивается.

Как показано на Рис. 2.10 а, за 5 итераций норма невязки убывает до минимально возможного значения. Оно обусловлено конечной точностью ма шинных чисел.

Мы используем арифметику двойной точности, то есть 8-ми байтовые чис ла, имеющие около 16 десятичных разрядов.

Разность между приближенным Vh и точным V решениями можно харак теризовать ее нормой в пространстве L1(). Разделив эту норму на объем области ||, получаем более наглядное среднее значение модуля функции Vh V :

(2.84) Vh V = |Vh V | d.

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

Тогда интеграл превращается в сумму разностей узловых значений с веса ми, равными объемам ячеек. Чтобы исключить приграничную погрешность, можно рассматривать Vh V только в некоторой подобласти, проведя расче ты в более широкой области. Например, при рассматриваемых здесь расчетах для n = 100 в области 600 км 600 км в подобласти 400 км 400км полу чаем значения нормы (2.84) 0.00759, 0.00200, 0.00052. Уменьшение примерно вчетверо при мельчении сетки вдвое демонстрирует близкую к квадратичной сходимость приближенных решений к точному. Отметим, что сама функция V имеет значения порядка единицы.

На Рис. 2.10 б показано жирными ломаными показано убывание этой нор мы погрешности потенциала Vh V в процессе итераций. Тонкими ломаны ми показано отличие Vh после данной итерации от того Vh, которое устано вится после 20 итераций. Как видим, сходимость процесса решения системы линейных алгебраических уравнений (2.29) в этой норме продолжается еще несколько итераций после того, как норма невязки уже достигла минимально го значения. После нескольких первых итераций погрешность приближенного решения Vh V уже не уменьшается, поскольку она обусловлена погрешно стью аппроксимации.

Рис. 2.10 б демонстрирует, что для самой грубой сетки итерации практи чески не уменьшают норму Vh V, она остается той же, что и была сразу после построения начального приближения. Это объясняется тем, что раз ность между точным и приближенным решениями V Vh 0.05, и разность между одномерным и точным решениями V V1D 0.08 отличаются всего в полтора раза, что мало заметно в масштабе Рис. 2.10 б. Однако малое изме нение нормы погрешности потенциала в процессе итераций не означает, что итерации не нужны. Достаточно обратить внимание на вертикальную ком поненту плотности тока jr. Начальное приближение, которое строится как одномерное решение, соответствует нулевому полиному Лежандра, и имеет практически постоянную по высоте jr = 1017 [А/м2 ]. В верхней атмосфере отличие такой jr от рассматриваемого точного решения, соответствующего сотому полиному Лежандра достигает четырех порядков, и почти устраня ется итерациями.

Для достижения практически той же точности всех представленных пара метров, что и после большого количества итераций, достаточно двух многосе точных итераций. На обычном персональном компьютере с процессором поко ления P4 и частотой 2 ГГц они занимают около 15 сек для сетки 16016032, которая имеет около миллиона узлов. Все массивы еще входят в оператив ную память и занимают 270 МБ. Гораздо больше времени, около 100 сек, занимает вычисление коэффициентов матрицы (2.29).

2.13. Параллельные вычисления.

Представленный в разделе 2.5. способ построения сетки позволяет дополни тельно разрезать любой криволинейный шестигранник. В описанном в раз деле 2.4. алгоритме основной объем вычислений требуется при проведении итераций метода верхней релаксации на основной и вспомогательных сетках.

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

Существует множество решений распараллеливания алгоритмов на мно гопроцессорной системе. Для научных и инженерных расчетов чаще всего применяются следующие модели параллельного программирования:

• Модель передачи сообщений. Каждый процесс обладает собственным локальным адресным пространством. Для синхронизации и обработки общих данных используется передача сообщений. Стандартом интер фейса передачи сообщений является MPI.

• Модель с общей памятью. Все процессы разделяют единое адресное про странство. Доступ к общим данным регулируется с помощью примити вов синхронизации. Стандартом для моделей с общей памятью стал OpenMP.

• Модель параллелизма по данным. В этой модели данные разделяются между узлами вычислительной системы, а последовательная програм ма их обработки преобразуется компилятором в программу либо в мо дели передачи сообщений, либо в модели с общей памятью. При этом вычисления распределяются по правилу собственных вычислений: каж дый процессор выполняет вычисления для данных, распределенных на него. Примером реализации этой модели является стандарты HPF1 и HPF2. На модели параллелизма по данным была также разработана отечественная система DVM.

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

В ИВМ СО РАН установлен кластер - аналог МВС100К, система с рас пределенной памятью. Эффективное использование такого типа многопро цессорных систем достигается при равномерной декомпозиции данных и вы числений по процессорам с минимальным количеством обменов информацией между процессами.

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

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

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

Алгоритм численного решения задачи, описанный в разделе 2.4., требу ет наибольшего количества вычислений внутри шестигранников. Мы имеем в виду шестигранники, на которые расчетная область разрезана, чтобы в каждом из них была правильная нумерация узлов сетки. Иногда мы их до полнительно разрезаем для достижения примерного равенство числа узлов в каждом, чтобы сделать более равномерным распределение вычислений по процессорам.

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

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

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

Более подробно работа параллельной части программного комплекса опи сана в Приложении A.

Таблица 1. Время расчетов в зависимости от количества используемых процессоров Количество процессоров 2 3 4 6 Время расчетов (сек) 545 155 116 80 Проверка эффективности параллельных вычислений проводилась для об ласти, описанной в разделе 2.5., и дополнительно разделенной по горизонта ли. Центральный шестигранник был разделен на четыре части, а каждый боковой – на две.

В каждом из получившихся 12 шестигранников была по строена сетка 256 256 16. Для решения всей задачи на одном компьюте ре потребовалось бы около 3 гигабайт оперативной памяти. Большой объ ем памяти требуется в нашем комплексе программ, поскольку он рассчитан на достаточно произвольные сетки, когда уравнения системы линейных ал гебраических уравнений вариационно-разностной схемы связывают значения потенциала в узле со значениями во всех 26 соседних узлах. Из-за симмет рии матрицы храним в памяти только 14 коэффициентов на узел. В табли це 1 приведено время расчетов в зависимости от количества используемых процессоров при проведении 10 многосеточных итераций. Для тестирования использовался кластер, установленный в ИВМ СО РАН, имеющий 48 процес соров с раздельной памятью по 1 гигабайту на процессор.

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

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

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

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

Это не было сделано, поскольку реализация кластера в ИВМ СО РАН на момент выполнения работы не позволяла использовать многоядерность на узлах кластера, так как каждое ядро представлялось отдельным узлом. Кро ме того, современные компиляторы эффективно векторизуют, оптимизируют и автоматически вставляют мультиядерное распараллеливание циклов. Хотя на данный момент директивы OpenMP позволяют получить более оптималь но распараллеленную программу, сам компилятор генерирует достаточно эф фективный код для процессора.

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

60o 40o 20o Рис. 2.11. Распределение Er на высоте 22 км – жирные линии. Между линиями уровня 13.5 мВ/м. Одномерное решение показано тонкими линиями, которые совпадают с эквипотенциалями в ионосфере с интервалом 5 кВ.

Соответствующие расчеты проведены для северного полушария на сет ке, представленной на Рис. 2.2, 2.3, но содержащей 513 513 17 узлов в центральном криволинейном шестиграннике и вдвое меньше - в остальных четырех. Поскольку основные поля сосредоточены в высоких и средних ши ротах, получаются те же результаты, если ограничиться одним центральным шестигранником, немного расширенным, чтобы покрыть интересующую нас область. Поскольку шаги сетки соответствуют выбранным в тестовых рас четах предыдущего параграфа, сходимость итераций и сходимость по шагу сетки не отличаются от тестовых. Сетка 513 513 17 содержит около 4. миллионов узлов. При независимых расчетах в каждой из четвертей этой области, расширенных на 100 км по горизонтали, программа умещается в области оперативной памяти 270 МБ, используемой имеющимся у нас ком пилятором Фортрана. Тогда реальное время работы компьютера примерно равно суммарному времени работы процессора, 3 минутам. Для всей области время работы процессора, как и размеры массивов, увеличивается вчетверо, а реальное время работы компьютера возрастает в зависимости от оператив ной памяти компьютера. При 1000 МБ – примерно до тех же 12 минут, а при 500 МБ – до 70 минут, то есть из-за обращений к внешней памяти получает ся шестикратное замедление. При увеличении числа узлов сетки аналогичное замедление происходит и на компьютере с большей памятью. Поэтому целе сообразно вести независимые расчеты в подобластях. Результаты расчетов, как и в описанных выше тестах, не меняются.

Существенным свойством исходного распределения потенциала на верх ней границе атмосферы является его гладкость. Если использовать исходное распределение потенциала на достаточно грубой сетке, и его интерполиро вать на расчетную сетку, результаты в верхней атмосфере существенно опре деляются выбранной интерполяцией. При линейной интерполяции имевши еся скачки производных на границах ячеек исходной сетки проявляются и у потенциала на расчетной мелкой сетке. Картина эквипотенциалей выгля дит довольно гладкой – тонкие линии на Рис. 2.11, но в разложении такой функции есть высокие гармоники. Как мы видели на Рис. 2.7, рост порядка полинома Лежандра дает быстрый рост Er, то есть вертикальной производ ной потенциала. Решение существенно меняется. Например, в горизонталь ных сечениях вертикальной компоненты электрического поля Er в верхней половине атмосферы утрированно вырисовываются все ребра исходной гру бой сетки, которая использовались для расчетов ионосферного поля в [62].

Аналогичный Рис. 2.11 вид сохраняется до тем больших высот, чем более гладким берется исходное распределение потенциала на верхней границе ат мосферы. Если провести сто итераций сглаживания исходного распределения потенциала, каждый раз заменяя значение в узлах на средние арифметиче ские значений в четырех соседних узлах, на высоте 60 км получаем поле, показанное на Рис. 2.12. При таком сглаживании увеличивается радиус кри визны эквипотенциалей в местах их наибольшего изгиба. Участки линий при этом смещаются на расстояния порядка расстояний между линиями. Под черкнем, что Рис. 2.11, 2.12 относятся к разным решениям, отличающимся сглаживанием исходных данных.

60o 40o 20o Рис. 2.12. Распределение Er на высоте 60 км – жирные линии. Между линиями уровня 69 мкВ/м. Одномерное решение показано тонкими линиями, которые совпадают с эквипотенциалями в ионосфере с интервалом 5 кВ.

Если рассматривать линии уровня Er на фиксированной большой высо те, где картина совершенно неудовлетворительна, со сглаживанием исходного потенциала исчезают мелкомасштабные погрешности, и картина становится гладкой, но сильно отличается от одномерной. Как видно на Рис. 2.12, поле имеет иное пространственное распределение, и его напряженность в обла стях экстремумов удвоена по сравнению с одномерным решением, показан ным тонкими линиями. Следует отметить, что грубая расчетная сетка тоже увеличивает долю верхней атмосферы, в которой сказываются погрешности задания потенциала на верхней границе.

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

2.15. Выводы.

Построен эффективный вычислительный алгоритм для решения глобальной задачи электропроводности атмосферы.

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

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

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

Глава Проникновение электрического поля из ионосферы до поверхности Земли 3.1. Введение.

При изучении атмосферных электрических полей основное внимание обычно уделяется грозовым явлениям. Глобальным результатом гроз является под держание разности потенциалов между поверхностью Земли и ионосферой около 300 кВ. Поскольку атмосфера является проводником, эта разность по тенциалов обеспечивает примерно постоянный электрический ток из ионо сферы на землю. Соответствующее электрическое поле вблизи поверхности Земли в среднем составляет 130 В/м. Заметим, что при высотном распре делении атмосферной проводимости соответствующем эмпирической модели [111], воспроизведенном на Рис. 1.1 а, такое поле получается при потенциале ионосферы около 350 кВ.

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

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

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

Интересен и обратный эффект – проникновение ионосферных электриче ских полей через атмосферу до земли. В частности, такие поля необходимо учитывать при объяснении наблюдаемых в приземном слое атмосферы вари аций электрического поля.

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

В настоящей главе мы используем модель ночной проводимости, пред ставленную на Рис. 1.1 а.

3.2. Направление распространения поля Поясним, в каком смысле мы говорим о проникновении электрического поля из ионосферы до поверхности Земли.

+ z, км а x, км 0 100 + + б в Ri Ra Ra Ra Ra Рис. 3.1. а - Распределение потенциала в вертикальном сечении атмосферы.

Эквивалентная электротехническая схема с генератором в ионосфере - б, или в магнитосфере - в. Стрелками показаны направления токов.

Положим атмосферу плоской и магнитное поле вертикальным. В условии (2.1) на верхней границе атмосферы зададим потенциал в виде (3.1) V |z=80 км = sin(kx), с периодом 200 км, и в соответствии с (2.2) полагаем потенциал нулевым на земле (3.2) V |z=0 = 0, Решим задачу электропроводности (1.17, 3.1, 3.2). Полученное распреде ление потенциала в вертикальном сечении атмосферы для одного периода по x показано на Рис. 3.1 а. Отметим, что это решение примерно соответствует гармонике n = 200, для которой точное решение было показано на Рис. 2.7.

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

2.7. На Рис. 3.1 а видно распределение потенциала не только по вертикали, как на Рис. 2.7, но и по горизонтали.

Стрелками на Рис. 3.1 показаны направления токов. Поскольку в ионо сферу положительный заряд вносится в области отрицательного потенциала и покидает ионосферу в области положительного потенциала, в ионосфере над зарядом совершается положительная работа, и значит там расположен источник энергии. На Рис. 3.1 б приведена эквивалентная электротехниче ская схема с генератором в ионосфере. Как и всегда при определении внут ренней структуры "черного ящика"по измерениям снаружи, при тех же элек трических полях и токах ниже ионосферы может быть и более сложная кон струкция генератора, например, показанная на 3.1 в. Именно она справедлива для рассматриваемого в настоящей главе случая, поскольку использованные нами ионосферные электрические поля получены в работе [62] как следствие токов из магнитосферы.

Поэтому справедливо говорить, что в построенной нами модели поля в атмосферу попадают из ионосферы, но генераторы расположены в магни тосфере. Если бы мы использовали как входные данные ионосферные элек трические поля, порождаемые движением ионосферного проводника в гео магнитном поле [31], то генератор был бы в ионосфере, и была бы эквива лентной электротехническая схема, приведенная на Рис. 3.1 б, но при расчете атмосферных полей мы решали бы ту же краевую задачу (1.17, 2.1, 2.2).

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

а б z, км x, км 0 100 0 Рис. 3.2. Распределения потенциала в вертикальном сечении атмосферы.

Гармонические по горизонтали решения с периодами x = 200 км – а, и 40 км – б. Между соседними эквипотенциалями V = 0.1. Тонкими линиями вблизи земли на фрагменте б дополнительно проведены эквипотенциали с шагом V = 0.005. Наклонные прямые показывают направление магнитного поля, которое отличается от вертикали на 300.

3.3. Влияние наклона магнитного поля Показанное на Рис. 3.1 а распределение потенциала получено при вертикаль ном магнитное поле. Сохраним все условия задачи, только зададим магнитное поле, отклоненное на 300 от вертикали. Положительный полупериод соответ ствующего решения показан на Рис. 3.2 а. Правый фрагмент этого рисунка демонстрирует поведение решения с горизонтальным периодом 40 км, кото рое близко к представленной на Рис. 2.7 гармонике n = 1000. Это решение относится к мелкомасштабным. Оно существенно затухает в верхней атмо сфере.

Как видно по этому рисунку, наклон магнитного поля приводит к сдвигу эквипотенциалей вдоль магнитных силовых линий на высоту около 65 км.

а б z, км x, км 0 100 0 Рис. 3.3. То же, что и Рис. 3.2, но при наклоне магнитного поля на 600.

Если высоту верхней границы, на которой задано ионосферное распре деление потенциала (3.1), увеличить до 120 км - уровня максимальной пе дерсеновской проводимости, то в масштабе рисунка ниже 80 км изменения не заметны. Следует только отметить, что в мелкомасштабном решении с периодом 40 км наряду с переносом вдоль магнитных силовых линий при снижении от 120 км до 80 км потенциал убывает примерно на 2%, а в круп номасштабном решении с периодом 200 км это уменьшение потенциала совсем мало – около 0.1%.

Аналогичный Рис. 3.3 построен при вдвое большем наклоне магнитного поля, 600 от вертикали. Как видно по этому рисунку, наклон магнитного поля приводит к сдвигу эквипотенциалей вдоль магнитных силовых линий на высоту около 70 км.

3.4. Результаты моделирования.

В работе [62] построены глобальные модели ионосферных электрических по лей, создаваемых продольными токами в авроральной зоне во время магни тосферных суббурь и в спокойных геомагнитных условиях. Здесь рассмотрим вклад этих полей в атмосферное электрическое поле.

Соответствующие пространственные распределения электрического поля получены в результате численного решения задачи электропроводности (1.6, 2.2, 2.1) при высотном распределении атмосферной проводимости, показан ном на Рис. 1.1 а. Поскольку сторонние токи в атмосфере мы не рассматри ваем, q = 0. На Рис. 3.4 приведены полученные распределения вертикальной компоненты на поверхности Земли. Их отличия от результатов одномерной модели не видны в масштабе рисунка. Это, в частности, означает совпадение приведенных линий с линиями уровня потенциала в ионосфере, построенны ми с шагом около 9 КВ.

Во время взрывной фазы интенсивной суббури, когда разность потенциа лов в ионосфере может достигать 120 кВ, на поверхности Земли получилось вертикальное поле от -32 В/м в вечернем секторе авроральной зоны до + В/м в утреннем секторе. В моделях [125] при разной геомагнитной активно сти разность потенциалов в ионосфере не превосходит 100 кВ, и соответству ющее вертикальное поле на поверхности Земли изменяется на 40 В/м. При меньшей геомагнитной активности получаются меньшие возмущения атмо сферного поля, однако даже в спокойных условиях они составляют ±10 % (а) (б) (в) Рис. 3.4. Распределения вертикальной компоненты электрического поля на поверхности Земли с интервалом 4 В/м. Направление на Солнце - вверх.

Пунктиром показаны отрицательные значения. а - спокойные условия, б взрывная фаза суббури, в - восстановительная фаза суббури.

среднего поля.

На Рис. 3.5 приведены распределения вертикальной компоненты электри ческого поля во время взрывной фазы суббури на высотах 35 и 50 км в срав нении с результатами одномерной модели. Как видим, различия существенно нарастают именно на этих высотах. Распределения на больших высотах мало похожи на получающиеся в одномерной модели.

400 Рис. 3.5. Распределения вертикальной компоненты электрического поля во время взрывной фазы суббури на высоте 35 км – (а), между линиями уровня 3 мВ/м. (б) – то же на высоте 50 км, интервал 1 мВ/м. Жирные линии – численное решение трехмерной задачи электропроводности, тонкие линии – одномерная модель.

3.5. Сопоставление результатов моделирования и измерений В работе [22] приведены результаты долговременных измерений вертикаль ной компоненты атмосферного электрического поля Ez на станции Восток в Антарктиде и выполнено сопоставление Ez с распределениями потенциала в ионосфере. При сопоставлении использовалась одномерная модель Ez (z), что для приземной атмосферы вполне оправдано, как это продемонстрировано на Рис. 2.9, на котором видно совпадение всех крупномасштабных решений с од номерным. Поэтому было достаточно знать лишь значение потенциала в той точке ионосферы, которая в данный момент находилась над станцией.

Для одномерной модели справедлива высотная зависимость потенциала (2.80), которая принимает вид z zup dz dz (3.3) V (z) = Vi, (z ) (z ) 0 если пренебречь кривизной земли и обозначить через Vi значение потенциала в ионосфере.

Для простоты интерпретации предположим экспоненциальный рост про водимости с высотой (3.4) (z) = 0 exp (z/h), где 0 и h – некоторые константы. Тогда из (3.3) с учетом большой толщины атмосферы, zup h получаем (3.5) V (z) = Vi (1 exp (z/h)), и вертикальное электрическое поле у поверхности Земли (3.6) Ez = Vi /h.

Для справедливости этого простого соотношения с погрешностью менее 10% достаточно, чтобы проводимость росла по закону (3.4) до z 2h и про должала расти выше этого слоя.

Распределения ионосферного потенциала в работе [22] задавались с помо щью двух моделей, Лукьяновой и др. [87] и Веймера [124]. Соответствующие им параметры мы отмечаем индексами L и W. Основные результаты, пред ставленные в работе [22] на Рис. 3, там же аппроксимированы линейными L W зависимостями, которые дают изменения поля Ez = 80 В/м и Ez = В/м при изменении значения потенциала V = 10 КВ. Отметим, что мы ис пользуем проведенные на этом рисунке прямые, а не их уравнения, которые приведены и на рисунке, и в тексте с очевидными опечатками.

В соответствии с (3.6) эти соотношения означают hL = 0.13 км и hW = 0. км. При таких малых величинах h потенциалу Vi = 300 КВ во всей ионосфере L W соответствуют поля у земли Ez = 2400 В/м и Ez = 1500 В/м, которые на порядок превосходят обычное поле хорошей погоды, которое принято считать равным 130 В/м. К сожалению, в статье [22] не приведено значение этого параметра в результатах их измерений, а из анализируемых данных среднее значение заранее вычтено.


Отметим, что в наших моделях, описанных в предыдущем разделе, Ez = 130 В/м получается, если потенциал ионосферы Vi = 350 КВ. В этом от ношении использованная нами модель атмосферной проводимости [111], вос произведенная на Рис. 1.1 а, эквивалентна экспоненциальной модели (3.4) с параметром h = 2.7 км.

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

При проведении таких исследований, как представленные в работе [22], есть несколько сложных этапов: измерение Ez, получение текущих значений ионосферного потенциала, а также способ обработки данных. Например, если обратиться к результатам измерений в течение конкретных четырех суток, представленных в [22] на рис. 2, и в качестве изменений поля Ez и потен циала V взять суточные вариации, то получаются значения h от 0.5 до км, то есть впятеро больше, чем дает проведенный выше анализ рис. 3, где представлено гораздо больше данных. Эти значения h гораздо ближе к тра диционным представлениям.

3.6. Влияние приземных неоднородностей проводимости Большое влияние на приземное электрическое поле оказывают неоднород ности атмосферной проводимости. В конкретные моменты времени проводи мость может существенно отличаться от использованных выше усредненных значений. По данным [75] характерно снижение проводимости около поверх ности Земли до 0.3 · 1014 См/м при увеличении концентрации аэрозольных частиц с радиусом 0.25 мкм до 1.5 · 108 м3 или повышение до 20 · 1014 См/м за счет повышения концентрации радона до 2 · 108 м3.

Характерной толщиной приземных слоев возмущения проводимости яв z, км 1015 1014 1013, См/м Рис. 3.6. Вертикальное распределение проводимости. Штриховая линия – эмпирическая модель [111]. Сплошные кривые – понижение за счет загрязнения атмосферы пылью или повышение за счет выброса радона.

ляется z0 = 250 м [75]. Мы взяли распределение возмущение проводимости по высоте в виде = 0 exp ((z/z0)2/2), (3.7) где 0 = 1014 См/м или 20 · 1014 См/м, чтобы примерно получились указанные в [75] значения при добавлении к невозмущенной величине 1.3 · 1014 См/м. Такие распределения показаны на Рис. 3.6. Штриховая линия соответствует эмпирической модели [111]. Горизонтальные размеры областей возмущения могут быть весьма различными, но мы для определенности взяли их равным вертикальному в одном направлении, x0 = z0, и втрое большим в другом, y0 = 3z0, с аналогичными (3.7) функциями. Для получения фонового поля -130 В/м задали разность потенциалов между Землей и ионосферой кВ.

В результате численного решения задачи электропроводности (1.6, 2.1, 2.3) получены пространственные распределения электрического поля. По скольку сторонние токи в атмосфере мы не рассматриваем, q = 0. На Рис. 3. приведены соответствующие распределения вертикальной компоненты на по верхности Земли и распределения электрического потенциала в центральном вертикальном сечении, y = 0. Видно увеличение от фонового поля -130 В/м до -246 В/м в центре облака пыли и уменьшение до -27 В/м в центре облака б а y, км г в z, км 1. 1.5 0 1. 0 1.5 x, км x, км Рис. 3.7. Электрическое поле при наличии в атмосфере облака пыли (а, в) или радона (б, г). а, б – распределения вертикальной компоненты электрического поля на поверхности Земли. Между линиями уровня В/м. в, г – распределения электрического потенциала в вертикальном сечении y = 0. Между линиями уровня 10 кВ.

радона. Как видим, неоднородности атмосферной проводимости существенно возмущают приземное электрическое поле.

Поскольку направление тока перпендикулярно эквипотенциалям, Рис. 3. (в) демонстрирует не только увеличение напряженности поля в области низ кой проводимости, но и обтекание этой области, а Рис. 3.7 (г) – уменьшение напряженности поля и концентрацию тока в области высокой проводимости.

Соответствующие плотности вертикального тока в центре облаков меньше фонового в 2.3 раза или больше в 3.5 раза. Отметим, что это делает неприме нимой одномерную модель, которая основана на независимости вертикальной плотности тока от высоты. Однако ее использование позволяет существенно уменьшить размер расчетной области. Так, результаты этого раздела могут быть получены в области |x| 2, |y| 3, z 2 км c погрешностью около 1%, если на всей границе задать одномерное распределение потенциала, показан ное кривой n = 0 на Рис. (2.10), с сохранением вертикальной компоненты фоновой напряженности электрического поля у земли.

3.7. Выводы.

Таким образом, проникновение электрического поля из ионосферы до по верхности Земли позволяет объяснить изменение вертикальной компоненты электрического поля вблизи земли примерно на 50 В/м, от 110 В/м до 160 В/м, при среднем наблюдаемом поле 130 В/м.

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

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

Измерения в Антарктиде показали на порядок больший вклад ионосфер ных разностей потенциала в вариации электрического поля около земли, чем получаются в нашей модели. Это несоответсвие можно было бы объяснить более резким возрастанием проводимости воздуха с высотой, чем в использу емой нами усредненной модели атмосферной проводимости [111]. Существует ли такой приземный слой в Антарктиде, могут сказать только натурные из мерения проводимости.

Гораздо сильнее приземное электрическое поле возмущается неоднород ностями атмосферной проводимости. В построенных моделях с типичными для атмосферы параметрами фоновое поле -130 В/м увеличивается почти вдвое в центре облака пыли или уменьшается впятеро в центре области с повышенным содержанием радона.

Исследована роль отличия направления магнитного поля от вертикаль ного. Показано, что переход от ионосферы с существенно тензорной прово димостью к атмосфере с практически скалярной проводимостью можно при близительно описать как перенос значений потенциала из Е-слоя ионосферы вдоль магнитных силовых линий на высоту около 70 км.

Глава Проникновение электрического поля от поверхности Земли в ионосферу 4.1. Предвестники сейсмической активности и их обнаружение со спутников. Обзор литературы.

Первые необычные возмущения вертикальной компоненты Ez электрическо го поля перед землетрясениям в области эпицентров на поверхности Земли наблюдались Кондо [84]. Флуктуации электрического потенциала с частотой десятки-сотни герц были обнаружены спутником на высоте 400 км [79]. В работе [72] приведены измерения квазистационарных возмущений электри ческого поля в ионосфере над областями, где предстоят или прошли земле трясения.

Кондо [84] предположил связь между выделениями радона перед земле трясениями и изменениям атмосферных параметров и электрического потен циала. Систематические исследования стали развиваться в 80-х годах [71].

Реакция атмосферы и ионосферы на интенсивные землетрясения была рас смотрена несколькими исследователями [78, 82, 36, 64, 112]. Наблюдались также возмущения, связанные с ядерными взрывами [113, 90], и связанные с распространением акустическо-гравитационных волн [110]. Детальный обзор этих и других исследований выполнен Парротом [99].

4.1.1. Модели предсказания сейсмической активности.

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

В своей работе Кондо [84] связывал это явление с увеличением освобожде ния радиоактивного газа - радона из земли перед и во время землетрясения.

В последующих работах в данном направлении King [81] и Pierce [101] при шли к выводу, что концентрация радона может быть увеличена искусственно или уменьшена плохой погодой.

Многие более поздние модели возникновения флюктуаций электрическо го поля основаны на возникновении микро-трещин за несколько дней перед землетрясением [93, 95]. Из-за механических сил в особенной части литосфе ры возникают микро-трещины, количество которых увеличиваются и в кон це концов приводит к землетрясению. Часть механической энергии, которая освобождается при землетрясении, превращается в электромагнитную энер гию [27, 97]. Получившееся электромагнитная эмиссия имеет те же часто ты, что и механические возмущения, и находится в низкочастотной области спектра (кГц и меньше). Поскольку некоторые части литосферы пропитаны водой, трещины также возникают при изменении давления воды в скважи нах, которые приводят к электрокинетическим эффектам [70]. Перемещения земной коры, например, из-за сейсмических волн, могут привести к индуктив ным эффектам, поскольку при этом перемещается проводник в геомагнитном поле.

Пьезо-электрическая теория возникновения электрических полей во вре мя землетрясений [66] была опровергнута измерениями проводимости земли [65], которая оказалась на несколько порядков меньше, чем необходимо, что бы работала пьезо-электрическая теория. Однако в настоящее время эта тео рия все же находит подтверждение в работе [67], в которой рассматривается возникновение токов в матрице минералов магматических и метаморфных пород под действием высокого давления, возникающего перед землетрясени ем.


Все эти эффекты, как и изменения в электромагнитном окружении око лоземной атмосферы из-за изменений ее химического состава могут привести к увеличению электромагнитного поля на поверхности Земли [69]. Несколь ко моделей предлагают суммировать результаты действия основных физи ческих механизмов и сопутствующих эффектов по всей цепочке, начина ющейся у поверхности Земли вплоть до ионосферы и даже магнитосферы [4, 106, 105, 116, 92]. Были найдены типичные изменения критической часто ты F2-слоя ионосферы foF2 [115], в полного электронного содержания [86], ионной температуры [114], локальной плотности ионов и электронов [100].

Было проведено множество работ по обработке и выявлению зависимостей между классом землетрясения, расстоянием до его эпицентра и изменением электрического поля [41].

Перечисленные механизмы образования аномалий скорее отражают слож ность быстро протекающих тектонических процессов в земной коре непосред ственно перед землетрясением, чем дают окончательный ответ о природе этих явлений.

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

Многие теоретические исследования этой проблемы базируются на рабо тах Дейнакаринтра и Парка [54]. Их подход был создан для изучения про никновения грозовых электрических полей в ионосферу и магнитосферу. Эта модель была развита и использована для улучшения описания физических явлений, происходящих во время разрядов молний [108, 98, 89, 120, 50, 80, 88, 121, 109]. Было выяснено, что проводимость атмосферы является ключе вым параметром в электрических процессах, происходящих между землей и ионосферой [76, 16, 77].

Второе применение модель Дейнакаринтра и Парка [54] нашла в иссле дованиях распространения предвестников сейсмической активности. В этом случае анализируется проникновение квазистационарного электрического по ля сейсмического источника в ионосферу [3, 17, 104]. Стоит также отметить модели, основанные на проникновение радиоактивных и заряженных газов в атмосферу перед землетрясением [101, 52, 117, 118]. Основным заключени ем перечисленных исследований является то, что электрическое поле может эффективно проникать в ионосферу и оказывать влияние на ионосферную плазму.

Однако существуют и статьи, в которых такая возможность опровергает ся [73, 56, 48]. В данной главе мы представим нашу модель и проанализируем указанное противоречие. Забегая вперед, отметим, что принципиальное раз личие результатов связано с корректностью учета ионосферного проводника.

Например, в моделях [104, 17] на верхней границе атмосферы было поставле но граничное условие, соответствующее нулевой проводимости ионосферы.

Проникновение поля более эффективно ночью, и величина напряженно сти поля существенно зависит от горизонтального размера источника, кото рый может быть областью готовящегося землетрясения [63, 107] заключили, что электрическое поле эффективно проникает в ионосферу, когда размеры исходной области порядка 100 200 км. Гримальский и др. [73] установили, что для существенного проникновения в ионосферу необходимо вертикальное приземное электрическое поле от 1 до 3 кВ/м. Такие значения обнаружены в результате измерений и обсуждаются в работе [91].

Широко известна еще одна модель проникновения электрического поля в ионосферу. Полученные на ее основе возмущения ионосферы представлены в статье [21], а две версии собственно модели усиления ионосферного элек трического поля – в статьях [42, 119]. От рассмотренных ниже трех моделей [56, 73, 17] эта модель отличается добавлением сторонних токов. Авторы по лагают, что существует сторонний вертикальный ток, представляющий собой вертикальное движение заряженных частиц аэрозоля, обусловленное конвек цией и турбулентной диффузией атмосферного воздуха. Однако такой сто ронний ток не может существовать как квазистационарный. В модели речь идет о переносе нескомпенсированных зарядах, но если в проводящую сре ду внести некоторый заряд, то как мы уже говорили в параграфе 1.1., он будет скомпенсирован за счет токов проводимости за время порядка 0/, и это время не превосходит 10 минут, поскольку проводимость обычно боль ше 2 · 1014См/м около поверхности Земли и быстро растет с высотой. Когда же положительные частицы будут окружены отрицательными частицами той же концентрации, перенос воздуха, содержащего все эти частицы, даст почти нулевой суммарный сторонний ток.

Необходимо сказать о существовании еще одного более простого объяс нения наблюдаемых в приземном слое атмосферы вариаций вертикального электрического поля – за счет изменения проводимости воздуха. Такие воз мущения поля подробно рассмотрены в последнем параграфе предыдущей главы, поскольку они обусловлены разностью потенциалов земля–ионосфера.

При этом не требуются какие-либо подземные электрические генераторы, ко торые являются причиной появления вертикального электрического поля в приземном слое атмосферы в моделях [73, 17] и в наших моделях, представ ленных в этой главе.

При наличии обычной разности потенциалов между Землей и ионосферой около 300 кВ изменение проводимости приводит к изменению токов и по лей, которые можно получить, решая уравнение электропроводности (1.17) с фиксированными значениями потенциала на верхней и нижней границах атмосферы. Если изменения проводимости не очень велики, и возмущением горизонтального электрического поля в ионосфере пренебречь, более общий подход не требуется. Как показано в работе [48], можно ограничиться одно мерной моделью, если горизонтальный масштаб явления не менее 100 км.

Именно такой упрощенный подход использован в работе [75]. В этой ста тье приведены зависимости проводимости от наличия радона и аэрозолей.

Так, десятикратное увеличение концентрации радона приводит к росту про водимости и уменьшает вертикальную компоненту электрического поля Ez в приземном слое атмосферы втрое, а десятикратное увеличение концентрации аэрозольных частиц с радиусом 0.25 мкм уменьшает проводимость и увели чивает это поле вчетверо. Соответствующие вариации Ez составляют + В/м и 400 В/м. Поскольку проводимость всегда остается положительной, этот эффект не может изменить знак Ez, но вариации Ez могут быть любо го знака и весьма велики. Например, большая положительная вариация Ez может быть, если сильно уменьшить загрязнение воздуха.

Однако выводы о существенном влиянии этого эффекта на ионосферу, сделанные на основе таких оценок в работе [75], представляются слишком смелыми, поскольку без обоснования подразумевается совпадение масштабов возмущения токовой скорости электронов в D слое ионосферы и скоро сти перемещения профиля ионизации D слоя, а эти скорости могут сильно различаться, так как ушедшие вверх электроны замещаются приносимыми снизу, и идут сложные процессы ионизации и рекомбинации. Следует также учесть, что эффект будет значительно слабее, если горизонтальный размер области выброса радона или аэрозоля меньше 100 км, поскольку задача элек тропроводности атмосферы перестает быть одномерной, и изменение полного тока в небольшой области вблизи Земли в нижней ионосфере перераспределя ется по гораздо более широкой области. При решении аналогичной теряющей одномерность задачи о проникновении токов от Земли в ионосферу в работе [48] показано, что при уменьшении до 10 км ширины областей положитель ного и отрицательного направления вертикального тока вблизи поверхности земли основная часть тока замыкается в нижней атмосфере, и доля этого то ка, замыкающаяся через ионосферу, снижается до менее чем одной десятой с существенным увеличением горизонтального размера областей втекания и вытекания тока. Соответственно, в десятки раз уменьшаются токовые ско рости электронов в D слое по сравнению со значениями из одномерных оценок.

Целью настоящей главы является изучение возможности проникновение приземного электрического поля через атмосферу в ионосферу. В этом случае именно приземное электрическое поле рассматривается как причина явления.

Основные результаты этих исследований опубликованы в работах [56, 12, 61].

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

4.2. Однослойная модель атмосферы В настоящем параграфе для наглядности сравнения с результатами ранних работ [17] рассмотрим максимально упрощенную модель, чтобы решения по лучались с использованием несложного математического аппарата.

Используем декартовы координаты, построенные в параграфе 1.4. Гео магнитное поле полагается вертикальным. Атмосфера рассматривается как горизонтальный слой между землей и некоторой плоскостью, расположен ной на высоте zup. Проводимость считается изотропной во всей атмосфере и зависящей только от высоты z. Высотное распределение аппроксимируем экспоненциальной функцией (4.1) (z) = exp (z/h), где константа равна значению у земли, и параметр h определяет скорость убывания с высотой. В соответствии с (Handbook of Geophysics, 1960) типич ные значения = 1013 С/м, h = 6 км для основной части атмосферы и = 3 · 1014С/м, h = 3км для области около Земли. Здесь мы используем z Ионосфера z = z up x=-a x=a x Земля Рис. 4.1. Однослойная модель атмосферы. Пунктирной кривой показано распределение вертикального электрического поля на поверхности Земли.

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

Выбор высоты zup, на которой заканчивается атмосфера и начинается ионосфера, не является тривиальным. Как видно по Рис. 1.1, чтобы про водимость во всей атмосфере можно было приближенно считать изотроп ной, желательно zup 70 км, а чтобы воспользоваться двумерной моделью ионосферного проводника, построенной в разделе 1.4. и основанной на суще ственной анизотропии проводимости P, надо взять zup 90 км. Мы выбрали zup = 80 км. Наши тестовые вычисления не показали значительных изменений результатов при выборе zup в пределах 80 90 км. Обоснование возможности такого условного разделения на атмосферу и ионосферу будет дано ниже в разделе 4.4. за счет использования более общей модели.

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

Геометрия модели показана на рис. 4.1. Аналогичный подход использовал ся в ранних работах [17, 104]. В данной двумерной модели уравнение элек тропроводности (1.17) принимает вид (4.2) (z) V (x, z) (z) V (x, z) = x x z z с нулевой правой частью, поскольку сторонние токи в атмосфере мы не рас сматриваем.

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

4.2.1. Нижнее граничное условие Как и в модели [104] полагаем, что по результатам измерений у поверх ности Земли задана вертикальная составляющая электрического поля. Тогда граничное условие для потенциала имеет вид (4.3) V (x, z) = E0(x), z z= где E0 (x) – заданная функция. Поскольку отсутствуют детальные простран ственные измерения, выберем гладкую функцию, отличную от нуля только в области |x| a (4.4) E0(x) = E0 (1 + cos(x/a))/2, |x| a, где a определяет размер области на поверхности Земли, и E0 — максимальное значение электрического поля, которое достигается в точке x = 0. Величина a может быть интерпретирована, как ширина области готовящегося земле трясения. Отрицательное вертикальное электрическое поле у земли означает увеличение электрического поля, существующего в спокойной атмосфере. В нашей модели были выбраны E0 = 100 В/м и a = 200 км. Эти значения соот ветствуют оценкам для процессов, предваряющих землетрясения умеренной магнитуды [91, 103].

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

Полагаем доли равными. Тогда исходные данные симметричны относитель но вертикали x = 0, и этим же свойством будет обладать решение.

4.2.2. Переход к ограниченной области От области, бесконечной направлении в x, можно перейти к конечной, не искажая решения вблизи источника, заданного условием (4.3), если удобным образом доопределить, как его ток замыкается на бесконечности. А именно, мы добавим (4.5) E0 = (1 + cos((x b)/a))/2, |x b| a, и продолжили эту функцию по всей оси x с периодом 2b. Величину параметра b a с помощью тестовых расчетов выберем достаточно большой, чтобы это доопределение не оказывало влияния на результат в интересующей нас области |x| 2a.

После этого условия задачи становятся периодическими, и следовательно этим свойством обладает решение:

(4.6) V (x + 2b, z) = V (x, z).

Из симметрии функции E0 следует также симметрия решения (4.7) V (b/2 x, z) = V (b/2 + x, z), V (x, z) = V (x, z), и поэтому потенциал равен нулю на вертикальных прямых x = ±b/ (4.8) V (±b/2, z) = 0.

Последнее условие позволяет ограничиться рассмотрением конечной об ласти |x| b/2 с вертикальными идеально проводящими стенками, достаточ но удаленными, чтобы не вносить погрешность в интересующей нас области |x| 2a. Такой прием, постановка некоторого искусственного граничного условия на некоторой достаточно удаленной границе, обычен для построения решений в бесконечных областях. В нашем случае он эквивалентен переходу к области |x| b с условием периодичности (4.6). С периодическими функ циями проще иметь дело, так как такие функции могут быть разложены в ряд Фурье, в то время как необходим интеграл Фурье (Korn and Korn, 1968) чтобы разложить функцию, которая определена на бесконечной оси и равна 0 снаружи конечного интервала |x| a как E0 (x) вида (4.4).

4.2.3. Верхнее граничное условие Мы рассматриваем ионосферу как слой zup z z, который описыва ем с помощью уравнения (1.16) для потенциала в ионосфере, полученного в приближении =.

Выражая плотность тока из атмосферы в ионосферу с помощью закона Ома (1.7), из (1.16) получаем V (x, y, zup) V (x, y, zup) (P ) (P )+ x x y y V (x, y, z) (4.9) + (zup) = 0.

z z=zup При таком выражении правой части уравнением (1.16) мы подразумеваем, что тока из ионосферы в магнитосферу нет V (4.10) (zup) = 0.

z z=z Поэтому Q в уравнении (1.16) состоит только из тока в ионосферу из ат мосферы, который заранее не известен, и поэтому в записи (4.9) перенесен в левую часть. Более общий случай, когда токи из ионосферы в магнито сферу существуют, например, из-за наличия сопряженной ионосферы, будет рассмотрен ниже. Отметим, что условие сопряжения решений в атмосфере и ионосфере (4.9) есть закон сохранения заряда, записанный с учетом локаль ного закона Ома для атмосферы и интегрального - для ионосферы.

С учетом независимости всех параметров от координаты y и скалярного характера проводимости (z) условие сопряжения (4.9) принимает вид V V (4.11) (P ) + (zup) = 0.

x x z z=zup z=zup Мы используем только постоянную интегральную проводимость P, так как горизонтальный масштаб интересующей нас области много меньше, чем характерный горизонтальный масштаб в ионосфере. Типичными являются значения P = 10 См в дневное время и в авроральной зоне и P = 0. См в ночное время. Примерно эти значения получается интегрированием ло кальной педерсеновской проводимости, представленной на Рис. 1.1 только по E слою ионосферы, то есть до высоты 150 км. Такое исключение прово димости F слоя ионосферы, как уже отмечалось в разделе 1.2., является традиционным при моделировании крупномасштабных электрических полей в ионосфере [31].

4.2.4. Решение краевой задачи Новая функция E1 (x) (4.5) равна E0 (x) (4.2) на интервале |x| b a и может быть представлена как:

E1(x) = fn cos(kn x), n= где kn = (2n 1)/b.

Члены с четными значениям 2n, так же как члены с sin(kx) отсутствуют, так как функция E1(x) - антисимметрична относительно x = b/2 и симмет рична относительно x = 0.

Коэффициенты Фурье равны:

2b fn = E1 (x) cos(knx)dx.

b Поскольку функция E1(x) симметрична относительно x = 0 и антисим метрична относительно x = b/2 антисимметрична относительно x = b/2, и базисные функции cos(knx) тоже обладают этими свойствами, можно вчет веро сократить интервал интегрирования b/ fn = E1(x) cos(knx)dx.

b Интервал интегрирования может быть уменьшен до 0 x a, так как E1(x) = 0 при a x b/2. Поскольку E1 (x) = E0(x) на интервале 0 x a может быть использована формула (4.2), и интегралы могут быть вычислены аналитически:

a 2 2 sin(kna) fn = E0(1 + cos(x/a)) cos(knx)dx = E0.

bkn((kna/)2 1) b Если kn a =, что означает 2n0 1 = b/a, то числитель и знаменатель равны нулю, и тогда fn = a/b. Этого может не случиться, если b/a - не целое.

Коэффициенты fn велики для n близких к n0 = (b/a1)/2 и уменьшаются как 1/n3 при n.

Общее решение уравнения (1.17) может быть найдено методом разделения переменных. Обычным образом получаем частные решения в виде произве дения функций, зависящей только от x, и зависящей только от z:

(4.12) V (x, z) = cos(kn x)(An exp n z + Bn exp n z), где 1 ) + k2, n = + ( 2h 2h 1 ) + k2.

n = ( 2h 2h Уравнения для коэффициентов An и Bn получаются из граничного усло вия (4.3) и условия (4.11). Обозначив отношение P k 2 + exp(n zup) n = exp ((n n )zup), P k 2 + exp(n zup ) получаем следующие выражения для коэффициентов An = fn /(n n n ), Bn = n An.

При найденных коэффициентах значение потенциала V (x, z) может быть вычислено в любой точке по формуле (4.12). В силу (1.5) компоненты элек трического поля могут быть получены почленным дифференцированием ря да (4.12) Ex (x, z) = V (x, z), x (4.13) Ez (x, z) = V (x, z), z а третья компонента Ey (x, z) = 0, так как потенциал не зависит от y.

4.2.5. Полученные результаты, См/м 1 9 Ez E, В/м Рис. 4.2. Высотные распределения вертикальной компоненты электрического поля и атмосферной проводимости. Пунктиром изображены распределения для однослойной модели над точкой x = 0. Сплошные кривые – для двухслойной модели над точкой x = amax.

Результаты расчетов электрического поля представлены пунктиром на Рис. 4.2, где показано высотное распределение вертикальной компоненты электрического поля на линии x = 0. Этот график почти не отличается от графика функции (4.14) Ez (0, z) = Ez (0, 0)(0)/(z), соответствующей не изменяющейся с высотой плотности вертикального элек трического тока.

E z, x, км В/м Рис. 4.3. Распределения вертикальной компоненты электрического поля на горизонтальных прямых, полученные в однослойной модели. Жирной кривой показано поле у земли, z = 0 км, тонкой – поле на высоте z = zup = 80 км, умноженное на exp(zup /h). Максимальные значения равны 100 В/м и 0.15 мВ/м, соответственно.



Pages:     | 1 || 3 | 4 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.