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

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

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


Pages:   || 2 | 3 | 4 |
-- [ Страница 1 ] --

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ

УЧРЕЖДЕНИЕ НАУКИ

ИНСТИТУТ ВЫЧИСЛИТЕЛЬНОГО МОДЕЛИРОВАНИЯ

СИБИРСКОГО ОТДЕЛЕНИЯ РОССИЙСКОЙ АКАДЕМИИ НАУК

На правах рукописи

ПОМОЗОВ

Егор Владимирович

УДК 519.632 + 551.594

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

КВАЗИСТАЦИОНАРНЫХ ЭЛЕКТРИЧЕСКИХ

ПОЛЕЙ В АТМОСФЕРЕ ЗЕМЛИ

Специальность 05.13.18 – Математическое моделирование, численные методы и комплексы программ ДИССЕРТАЦИЯ на соискание ученой степени кандидата физико-математических наук

Научный руководитель:

д.ф.-м.н. Денисенко В.В.

Красноярск – 2013 Оглавление Глава 1. Квазистационарная модель электрических полей и токов в атмо сфере и ионосфере Земли 1.1. Основные уравнения........................... 1.2. Проводимость............................... 1.3. Двумерная модель ионосферного проводника............. 1.4. Запись уравнений в декартовых координатах............. Глава 2. Численный метод решения уравнения электропроводности со ска лярной проводимостью 2.1. Введение................................. 2.2. Формулировка задачи.......................... 2.3. Энергетический метод.......................... 2.4. Вариационно - разностная схема..................... 2.5. Построение сетки............................. 2.6. Многосеточный метод........................... 2.7. Сетка с кубическими ячейками..................... 2.8. Разностная схема............................. 2.9. Дискретная модель............................ 2.10. Сходимость к точным решениям..................... 2.11. Точные решения для сферически симметричной атмосферы..... 2.12. Фактическая сходимость к точным решениям............. 2.13. Параллельные вычисления........................ 2.14. Рекомендации по организации расчета атмосферного электриче ского поля................................. 2.15. Выводы................................... Глава 3. Проникновение электрического поля из ионосферы до поверхно сти Земли 3.1. Введение.

................................. 3.2. Направление распространения поля.................. 3.3. Влияние наклона магнитного поля................... 3.4. Результаты моделирования........................ 3.5. Сопоставление результатов моделирования и измерений....... 3.6. Влияние приземных неоднородностей проводимости......... 3.7. Выводы................................... Глава 4. Проникновение электрического поля от поверхности Земли в ионосферу 4.1. Предвестники сейсмической активности и их обнаружение со спут ников. Обзор литературы......................... 4.1.1. Модели предсказания сейсмической активности........... 4.1.2. Известные модели проникновения квазистационарного электри ческого поля от поверхности Земли в ионосферу.......... 4.2. Однослойная модель атмосферы.................... 4.2.1. Нижнее граничное условие...................... 4.2.2. Переход к ограниченной области................... 4.2.3. Верхнее граничное условие...................... 4.2.4. Решение краевой задачи....................... 4.2.5. Полученные результаты........................ 4.3. Двухслойная модель атмосферы.................... 4.3.1. Биполярный источник......................... 4.3.2. Проводимость............................. 4.3.3. Условия сшивки............................ 4.3.4. Переход к ограниченной области................... 4.3.5. Решение краевой задачи....................... 4.3.6. Результаты моделирования...................... 4.4. Возможность использования двумерной модели ионосферного про водника.................................. 4.5. Анализ погрешностей моделирования атмосферного электрическо го поля, возникающих при использовании упрощенных моделей ионосферного проводника........................ 4.6. Выводы.................................. Заключение Литература Приложение A. Программный комплекс A.1. Структура программного комплекса.................. A.2. Алгоритм работы однопоточной single-программы.......... A.2.1. Single-программа............................ A.2.2. Многосеточный метод......................... A.2.3. Метод последовательной верхней релаксации............ A.3. Алгоритм работы параллельной программы............. A.3.1. Подготовка данных для параллельной программы......... A.3.2. Multiproc-программа.......................... A.3.3. Обработка полученных данных.................... Приложение B. Построение сетки в расчетной области B.1. Элементы сетки.............................. B.2. Типы прилегания............................. B.3. Переменные................................ B.4. Массивы.................................. B.5. Нумерация узлов сетки......................... B.6. Сеточные функции............................ B.7. Алгоритм построения сетки....................... Приложение C. Фактическая сходимость численного метода C.1. Задача в параллелепипеде........................ C.2. Задача в призме с шестью боковыми гранями............ C.3. Задача в шаровом слое......................... Введение Объектами исследования являются квазистационарные электрические поля, возникающие в атмосфере Земли за счет генераторов, расположенных в магнитосфере и в литосфере.

Актуальность исследований.

Математическому моделированию крупномасштабных электрических по лей в атмосфере Земли посвящены многочисленные работы [3, 96, 51, 24].

Наряду с анализом полей, создаваемых грозовыми облаками [54, 121], ин тенсивно исследуются процессы, в которых атмосфера является связующим звеном между ионосферой и литосферой [17, 23, 119, 45]. Прикладная направ ленность моделирования этих процессов обусловлена желанием использовать космические средства для обнаружения предвестников землетрясений. Име ются экспериментальные данные об изменениях электрического поля в при земной атмосфере накануне землетрясений. Поскольку покрыть Землю до статочно плотной сетью наземных датчиков проблематично, возник вопрос, можно ли судить об этих полях на основе спутниковых измерений. Возмуще ния ионосферы накануне землетрясений тоже наблюдаются. Так, флуктуации электрического потенциала с частотой десятки-сотни герц были обнаружены спутником на высоте 400 км [79]. В работе [72] приведены измерения квази стационарных возмущений электрического поля в ионосфере над областями, где предстоят или прошли землетрясения.

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

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

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

Для достижения поставленной цели выделены следующие задачи:

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

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

Положения, выносимые на защиту:

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

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

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

Данные положения соответствуют трём пунктам паспорта специальности 05.13.18 — "Математическое моделирование, численные методы и комплексы программ"по физико-математическим наукам:

1. Разработка новых математических методов моделирования объектов и явлений.

4. Реализация эффективных численных методов и алгоритмов в виде комплексов проблемно-ориентированных программ для проведения вычис лительного эксперимента.

7. Разработка новых математических методов и алгоритмов интерпрета ции натурного эксперимента на основе его математической модели.

Личный вклад автора.

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

Автор построил математические модели крупномасштабных электриче ских полей, существующих в атмосфере Земли.

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

Научная новизна работы определяется тем, что:

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

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

проверкой и обоснованием точности использованных численных методов;

сравнением полученных результатов с известными в научной литературе теоретическими и экспериментальными результатами других авторов;

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

Теоретическая и практическая значимость.

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

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

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

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

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

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

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

По теме диссертации опубликовано 9 статей, из них 3 работы – в жур налах, рекомендуемых ВАК для защиты кандидатских диссертаций.

1. Денисенко В. В. Расчет глобальных электрических полей в земной ат мосфере. / В. В. Денисенко, Е. В. Помозов. – Вычислительные технологии.

Т. 15. N 5. 2010. – С. 34-50.

2. Denisenko V. V. Ionospheric conductivity eects on electrostatic eld penetration into the ionosphere. / V. V. Denisenko, M. Y. Boudjada, M. Horn, E. V. Pomozov, H. K. Biernat, K. Schwingenschuh, H. Lammer, G. Prattes, and E. Cristea. – Nat. Hazards Earth Syst. Sci. Vol. 8. 2008. – P. 1009–1017.

3. Denisenko V. V. On electric eld penetration from ground into the ionosphere. / V.V. Denisenko, M. Ampferer, E.V. Pomozov, A.V. Kitaev, W. Hausleitner, G. Stangl, H.K. Biernat– Journal of Atmospheric and Solar Terrestrial Physics. Vol. 102. 2013. P. 341-353.

В других изданиях:

4. Денисенко В. В. Расчет атмосферных электрических полей, проника ющих из ионосферы. / В. В. Денисенко, В. В. Бычков, Е. В. Помозов. – Солнечно-земная физика. Вып. 12. Т. 2. Иркутск: Изд-во ИСЗФ РАН. 2008.

– С. 281-283.

5. Денисенко В. В. Математическое моделирование проникновения элек трических полей из ионосферы в атмосферу. / В. В. Денисенко, В. В. Бычков, Е. В. Помозов. – Межгеосферные взаимодействия (Москва 26-27 сентября 2011 г.): материалы семинара-совещания / Ин-т динамики геосфер РАН. М.:

ГЕОС. 2011. – С. 89-96.

6. Denisenko V. V. Calculation of Atmospheric Electric Fields Penetrating from the Ionosphere. / V. V. Denisenko, V. V. Bychkov, E. V. Pomozov. – Geomagnetism and Aeronomy. 2009. Vol. 49. No. 8. – P. 1275-1277.

7. Денисенко В. В. Проникновение электрического поля из приземного слоя атмосферы в ионосферу. / В. В. Денисенко, Е. В. Помозов. Солнечно земная физика. – Вып. 16. Иркутск: Изд-во ИСЗФ РАН. 2010. – С. 70-75.

8. Denisenko V. V. Mathematical simulation of the large scale electric elds and currents in the Earth’s atmosphere. / V. V. Denisenko, E. V. Pomozov.

– Proceedings of the 8-th International Conference on Problems of Geocosmos, St.-Petersburg, September 20-24, 2010. – P. 392-397.

9. Denisenko V. V. Penetration of Electric Field from the Surface Layer to the Ionosphere. / V. V. Denisenko, E. V. Pomozov. – Geomagnetism and Aeronomy.

2011. Vol. 51. No. 7. – P. 866-872.

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

В трудах конференций опубликованы указанные выше статьи [5, 8], а также:

1. Помозов Е.В. Программный комплекс для решения трехмерных эллип тических уравнений многосеточным методом с использованием параллель ных вычислений. / Е. В. Помозов. – Материалы VII Всероссиийской конф.

молодых ученых по мат. моделированию и информационным технологиям.

ИВТ СО РАН, Новосибирск, 2006. – С. 25.

2. Помозов Е. В. Программный комплекс для нахождения электрическог поля в атмосфере многосеточным методом с использованием параллельных вычислений. / Е. В. Помозов. – Материалы конф. молодых ученых КНЦ СО РАН, Красноярск: КНЦ СО РАН, 2008. – С. 42-43.

В целом диссертация была доложена на семинаре отдела Вычислительной математики Института вычислительного моделирования СО РАН.

Структура и объем работы.

Диссертационная работа состоит из введения, четырех глав, заключения, списка литературы и трех приложений. Общий объем диссертации – 201 стра ница. Список литературы включает 125 наименований.

Содержание работы Во введении кратко описан объект исследования, сформулирована цель диссертационной работы, а также выделены ее задачи. Обоснована актуаль ность работы, теоретическая и практическая значимость, а также достовер ность полученных результатов. Перечислены основные публикации по теме диссертации и дана общая характеристика диссертации.

Поскольку в работе рассмотрены три существенно различающиеся задачи:

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

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

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

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

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

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

Существует ли такой специфический приземный слой в Антарктиде, могут сказать только натурные измерения проводимости.

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

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

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

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

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

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

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

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

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

Глава Квазистационарная модель электрических полей и токов в атмосфере и ионосфере Земли 1.1. Основные уравнения Поскольку мы интересуемся квазистационарными процессами в проводящей среде, основные уравнения для атмосферного электрического поля — это за кон Фарадея, закон сохранения заряда и закон Ома:

rot E = 0 (1.1) div j = q (1.2) (1.3) j = E, где E — напряженность электрического поля, j — плотность тока, — тензор проводимости. Заданная плотность источника тока q может быть отличной от нуля, если существуют сторонние токи. В наших моделях электрических полей мы не рассматриваем сторонние токи, но созданные методы решения задач пригодны и при их наличии. При заданной плотности сторонних токов jext (1.4) q = div jext.

Стационарная модель применима, если характерное время рассматрива емого процесса намного больше времени релаксации заряда = 0/. По скольку проводимость минимальна в приземном слое атмосферы, и См/м [94], в рамках такой модели можно рассматривать процессы, длящиеся намного более четверти часа.

Уравнение (1.1) позволяет ввести электрический потенциал, такой что (1.5) E = grad V.

Тогда система уравнений (1.1-1.3) сводится к уравнению электропровод ности div ( grad V ) = q, (1.6) которое превращается в уравнение Пуассона, если проводимость изотропна и однородна.

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

Закон Ома, который связывает электрическое поле E и плотность тока j, удобно записать, представив вектор E в виде суммы продольной компоненты E, параллельной вектору магнитного поля, и нормальной компоненты E.

Аналогичное представление используем и для j. Тогда j = E (1.7) j = P E H [E B] /B, где величины P, H, называются соответственно проводимостями Педер сена, Холла и продольной [1].

Для вычисления значений этих компонент тензора проводимости в ионо сфере выше 90 км мы используем модель, основанную на эмпирических мо делях пространственно-временных распределений ионосферных параметров IRI, MSISE, IGRF, представленную в работе [55]. В этой модели проводи мость рассчитывается до высоты 2000 км, однако мы используем распреде ления только до верхней границы F-слоя на высоте z = 500 км, поскольку выше лежащие слои добавляют менее 1 % к значениям интересующих нас параметров.

В большей части атмосферы Земли, ниже 60 км, тензор проводимости изотропен, так что H = 0, P =, и поэтому проводимость можно описы вать скалярным параметром, P = =. Его высотное распределение мы задаем в соответствии с эмпирической моделью [111].

В области высот от 60 до 90 км для P (z), (z) мы используем интер поляционные формулы, позволяющие гладко сшить эти модели. Поскольку проводимость Холла равна нулю ниже 60 км, в области интерполяции мы используем выражение (1.8) H (z) = P (z) (z) P (z), типичное для плазмы с одним преобладающим видом носителей заряда. Па раметр Холла, равный отношению H /P, в нижней ионосфере примерно ра вен отношению гирочастоты электронов к частоте столкновений электронов с нейтральными молекулами. Он принимает значения порядка 10 на высоте 90 км, и быстро убывает с уменьшением высоты. В нашей аппроксимации он обращается в нуль ниже 60 км, что соответствует изотропии тензора прово димости.

Типичные для средних широт высотные распределения компонент тензо ра проводимости, построенные нами на основе перечисленных выше эмпири ческих моделей, показаны на Рис. 1.1. Приведены дневные и ночные профи ли. Штриховыми линиями показаны эффективные педерсеновские проводи мости, которые позволяют учесть изменение плотности тока из-за ускорения ионосферной проводящей среды под действием силы Ампера в течение часа.

Детальное объяснение дано в работе [55]. Отметим только, что за длительное время среда может приобрести скорость дрейфа в скрещенных электриче ском и магнитном полях, и тогда напряженность электрического поля в си стеме покоя среды станет нулевой, и следовательно токи проводимости тоже станут нулевыми. Традиционным является упрощенный учет этого явления, состоящий в исключении проводимости ионосферного F слоя [31].

z, км Ночь День P H P H 1015 1010 105 1015 1010, См/м, См/м Рис. 1.1. Высотные распределения компонент тензора проводимости, построенные нами для средних широт. Ночные – слева, дневные – справа.

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

Следует отметить существование более сложных физических явлений в приземном слое атмосферы, которые не описываются в рамках используе мой простейшей физической модели. По-видимому, важнейшими из них яв ляются электродный эффект, заключающийся в дисбалансе положительных и отрицательных ионов вблизи поверхности Земли [83], а также формиро вание мелкомасштабных нестационарных электрических структур в этой об ласти [24, 49]. Мы также не рассматриваем сторонние токи, которые могут существовать за счет различия средних скоростей ионов разных знаков, обу словленного не электрическим полем, а, например, диффузией, как в модели [119].

1.3. Двумерная модель ионосферного проводника Отделим основную часть ионосферы от нижней ионосферы некоторой по верхностью и в основной части ионосферы будем использовать приближение =. Как видно на Рис. 1.1, выше 100 км на несколько порядков превосходит остальные компоненты тензора проводимости. Это приближе ние не справедливо в нижней ионосфере, и минимально возможная высота этой условной границы, при которой еще не вносится заметная погрешность, выбирается с помощью тестовых расчетов.

Математический аппарат, необходимый для перехода к двумерной модели в приближении бесконечной продольной проводимости, изложен в работе [5], и более простая версия для потенциального магнитного поля – в [55]. Опишем лишь основные приемы такого перехода.

Из = и конечности плотности тока следует равенство нулю про дольной компоненты электрического поля E = 0, и следовательно эквипо тенциальность магнитных силовых линий. Из равенства нулю продольной компоненты электрического поля и из определения потенциала (1.5) следует, что достаточно определить E на одной поверхности, пересекающей все си ловые линии магнитного поля. Такая поверхность названа базовой в работе [5]. Тогда во всех точках пространства оно может быть построено переносом вдоль силовых линий, то-есть для E модель становится двумерной.

В такой модели важную роль имеет интегральный ток между силовыми линиями, поскольку перераспределение заряда вдоль отдельной магнитной силовой линии происходит за счет бесконечно малого электрического поля E. Интегральный ток J определяется как интеграл от j вдоль лежащего в проводящем слое ионосферы отрезка силовой линии:

(1.9) J= j dl.

Так же можно проинтегрировать локальный закона Ома (1.7). При этом нормальные к магнитному полю компоненты напряженности электрического в каждой точке магнитной силовой линии с помощью несложных геометри ческих соображений выражаются через их значения в точке, лежащей на базовой поверхности, Eb Получаем двумерный закон Ома P H J= Eb (1.10) H P с педерсеновской P и холловской H интегральными проводимостями.

Если рассматриваемая область достаточно мала, чтобы магнитные си ловые линии рассматривать как параллельные, поле E просто неизменно вдоль каждой магнитной силовой линии E = Eb. Тогда его можно вынести из-под интеграла, и для интегральных проводимостей получаются простые выражения (1.11) P = P dl, H = H dl.

Закон сохранения заряда, проинтегрированный вдоль магнитной силовой трубки, означает, что дивергенция интегрального тока J компенсируется то ками через концы трубки, сумму которых обозначим через Q:

Div J = Q. (1.12) При заданной правой части это уравнение с использованием закона Ома (1.10) и выражением Eb через электрический потенциал V (1.5) становится двумерным уравнением электропроводности на базовой поверхности.

Такого рода модели обычно используются для расчета глобального элек трического поля в ионосфере. При этом правая часть (1.12) Q может содер жать заданные над ионосферой продольные токи из магнитосферы. Именно в результате таких расчетов в работе [62] получены распределения электри ческого потенциала в ионосфере, которые мы используем в Главе 3.

1.4. Запись уравнений в декартовых координатах В последней Главе 4 рассмотрены локальные явления, для которых можно не учитывать кривизну Земли. Если при этом магнитное поле вертикально, двумерная модель ионосферного проводника приобретает простой вид.

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

Определение интегрального тока J (1.9), который в рассматриваемом слу чае имеет x, y компоненты, принимает вид z z (1.13) Jx = jx dz, Jy = jy dz, zup zup двумерный закон Ома (1.10) – вид:

J H E x = P x, (1.14) Jy H P Ey а педерсеновская и холловская интегральные проводимости (1.11) получают ся интегрированием по высоте:

z z (1.15) P = P dz, H = H dz.

zup zup Уравнение электропроводности (1.12) для потенциала в ионосфере при нимает вид V V (1.16) P P = Q.

x x y y Отметим, что в рамках такой модели ионосфера эквивалентна тонкой про водящей пленке с двумерным законом Ома (1.14).

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

Уравнение электропроводности (1.6) для потенциала в атмосфере прини мает вид V V V (1.17) = q.

x x y y z z Глава Численный метод решения уравнения электропроводности со скалярной проводимостью 2.1. Введение Численное решение эллиптических краевых задач является одним из наибо лее развитых разделов вычислительной математики как в отношении реше ния прикладных задач, так и в теоретическом обосновании методов.

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

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

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

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

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

Однако для эллиптических краевых задач достаточно общего вида наи более эффективными представляются методы, использующие вариацион ные формулировки, например, принцип Дирихле для уравнения Пуассо на. Такие методы имеют много названий: методы конечных элементов [33], вариационно-разностные методы [32], проекционно-сеточные методы [26].

Во многих прикладных задачах вариационный принцип состоит в мини мизации значения функционала, совпадающего с выражением для энергии, и тогда методы называются энергетическими [28, 30]. В таких случаях пере ход к дискретной задаче упрощается по сравнению с общими вариационно разностными методами, поскольку из-за симметрии можно использовать ме тод Ритца вместо метода Бубнова-Галеркина. Существенным преимуществом энергетических методов является также симметрия и положительная опре деленность матриц систем линейных алгебраических уравнений, которые приходится решать, чтобы получить узловые значения аппроксимирующих функций.

Именно в рамках энергетического метода создан численный метод реше ния задачи электропроводности, представленный в настоящей главе.

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

Мы используем многосеточный метод Р.П. Федоренко [44], эффективность которого обусловлена как уменьшением числа арифметических операций при укрупнении ячеек сетки, так и с сужением диапазона собственных значений эллиптических операторов при укрупнении пространственного масштаба со ответствующих собственных функций. Многосеточный метод позволяет сде лать необходимое число итераций практически независимым от числа узлов сетки N, хотя теоретические оценки дают небольшое увеличение числа ите раций порядка ln N. Использование энергетического метода позволяет легко осуществлять переходы с мелких сеток на крупные с сохранением симметрии и положительной определенности матриц коэффициентов систем линейных алгебраических уравнений. Сочетание метода конечных элементов с много сеточным методом рассмотрено в монографии [46].

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

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

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

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

Стационарные электрические поля в атмосфере удовлетворяют уравне нию электропроводности (1.6), которое в декартовых координатах при изо тропной проводимости имеет вид (1.17).

Типичное высотные распределения компонент тензора проводимости, описанные нами в разделе 1.2, показаны на Рис. 1.1. В атмосфере ниже 60 км проводимость изотропна, поэтому - скаляр. В настоящей главе для демон страции возможностей созданного численного метода мы используем модель атмосферной проводимости, приведенную в обзоре [96].

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

Оба указанных упрощения не справедливы в слое, который является пе реходным между атмосферой и ионосферой, и обычно находится на высотах 7090 км. В этом слое нижней ионосферы проводимость с ростом высоты пе рестает быть изотропной, как в атмосфере, но проводимость в направлении B еще не становится много большей остальных компонент тензора прово димости, как в основных слоях ионосферы. Если исключить из рассмотре ния экваториальную область, где геомагнитное поле почти горизонтально, каждая магнитная силовая линия при прохождении рассматриваемого слоя, имеющего высоту около 20 км, смещается по горизонтали не более чем на несколько десятков километров. Поэтому, если экстраполировать эквипотен циальность магнитных силовых линий в этот слой, распределение ионосфер ного электрического потенциала, полученное на высоте 90 км, при переходе на высоту 70 км перераспределится по горизонтали на эти десятки километ ров. Это смещение можно рассматривать как оценку погрешности, вносимой тем, что в нашей модели мы заменяем нижнюю ионосферу поверхностью, на которой проводимость скачком меняется от изотропной атмосферной до существенно анизотропной ионосферной. Ниже этого скачка скалярную про водимость мы полагаем равной проводимости в вертикальном направлении, то есть при вертикальном геомагнитном поле, а не P, поскольку про водимость в горизонтальном направлении шунтируется хорошо проводящей ионосферой, и поэтому ее искажение не столь существенно. Разумеется, мы можем применять такую модель только для полей с горизонтальным масшта бом не менее сотни километров. Проведенные тестовые расчеты показывают, что выбор такой поверхности на высотах 80 90 км мало меняет получаю щееся атмосферное электрическое поле. Аккуратный учет этого переходного слоя может быть сделан в рамках трехмерной модели электропроводности с гиротропным тензором проводимости. Соответствующая модель проводимо сти приведена в [55], энергетический метод решения такой задачи, имеющей несамосопряженный оператор, предложен в [7]. В такой, более адекватной модели, целесообразно отделить атмосферу, добавив условия сшивки на но вой границе. Для атмосферы получится модель, используемая в настоящей работе. В силу изотропии проводимости методы решения для собственно ат мосферы существенно упрощаются.

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

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

(2.1) V = V0 (, ), r=Rup где r,, – сферические координаты, Rup – радиус сферы, которая исполь зуется в качестве верхней границы атмосферы.


Далее мы перейдем к задаче с однородным граничным условием (2.1), решение которой будем обозначать как V без тильды.

Нижней границей атмосферы является поверхность Земли. Характерная проводимость приземного воздуха порядка 1014 См/м, что на много поряд ков меньше проводимости морской воды, 3 См/м, влажной земли, 102 См/м, и даже таких пород, как мрамор, 108 См/м. Поэтому поверхность Земли обычно рассматривается как идеальный проводник, что соответствует гра ничному условию (2.2) V |r=RE = 0, где RE – радиус Земли.

Это условие записано в предположении сферичности Земли. Если учиты вать рельеф с заданной с помощью функции h(, ) высотой поверхности над уровнем моря, то условие Дирихле (2.2) принимает вид (2.3) V |r=RE +h(,) = 0.

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

n r=RE +h(,) Это условие может быть и неоднородным, с заданной нормальной к грани це компонентой электрического поля, как в моделях, рассматриваемых ниже в главе (4).

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

В отсутствии идеальных проводников и изоляторов внутри области про водимость удовлетворяет ограничениям (2.5) 0 1 2, где 1, 2 – константы.

В соответствии с [30] для разрешимости задачи Дирихле функция V0(, ) в граничном условии (2.1) должна быть такой, что в области существует функция (r,, ), которая на границах принимает заданные значения (2.6) |r=Rup = V0 (, ) (2.7) |r=RE +h(,) = и имеет конечную джоулеву диссипацию, то есть (grad )2d. (2.8) Для функций V0 (, ) с ограниченными первыми производными в каче стве может быть использована просто линейная интерполяция по радиусу от RE + h(, ) до Rup соответствующих значений, нуля и V0 (, ).

Наличие такой функции позволяет свести исследование неоднородной задачи Дирихле к задаче для функции V с нулевыми граничными усло виями. При численном решении задачи выделение некоторой функции мо жет иметь смысл начального приближения. Мы иногда используем функцию (r,, ), полученную как решение одномерной по r задачи, когда для каждо го луча, то есть при фиксированных,, строится сферически симметричное решение.

В предположениях (2.5, 2.8) для разности (2.9) V = V требуется решить однородную задачу Дирихле для уравнения электропро водности, получающегося из (1.6) :

div( grad V ) = q div( grad ) (2.10) (2.11) V |r=Rup = (2.12) V |r=RE +h(,) = 0.

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

2.3. Энергетический метод.

Для задачи Дирихле (2.10-2.12) справедлив принцип минимума функционала энергии [28], в соответствии с которым обобщенное решение краевой задачи может быть получено в результате минимизации функционала энергии (grad V )2 d 2 grad V · grad d, (2.13) W (V ) = qV d где · означает скалярное произведение векторов, интегрирование проводится по всей области, d – элемент объема.

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

Квадратичная часть W (V ) соответствует энергетическому скалярному произведению (2.14) [V, U ] = grad V · grad U d, где V, U – произвольные функции, удовлетворяющих условиям (2.11, 2.12).

Как доказано в [30] эта симметричная билинейная форма на рассматри ваемом множестве функций положительно определена:

(grad V )2d V 2 d, (2.15) c где c – некоторая константа.

Минимум функционала энергии существует, единственен и дает обобщен ное решение задачи (2.10-2.12), которое при дополнительном предположении гладкости является классическим решением.

(1) Энергетическая норма эквивалентна норме пространства W2 (). Поэто му обобщенное решение можно аппроксимировать кусочно-линейными функ циями [32].

2.4. Вариационно - разностная схема.

Для использования кусочно-линейных функций необходимо предварительно разбить область на тетраэдры. Считаем, что такое разбиение выполнено с со блюдением стандартных ограничений для нерегулярной сетки [32], одним из которых является требование, чтобы длинны всех ребер лежали в интервале l1h - l2h, где l1 и l2 - положительные константы, а параметр h называется шагом сетки.

Обозначим через I число внутренних узлов сетки. Через V h (x, y, z) обо значим искомую кусочно-линейную функцию, аппроксимирующую решение краевой задачи (2.10-2.12). В монографии [32] доказана в двумерном случае эквивалентность L2 нормы кусочно-линейной функции V h (x, y) и евклидо вой нормы набора ее узловых значений Vi (неравенства (5) из §2 главы 2).

Аналогичные неравенства для трехмерной функции V h (x, y, z) доказаны в монографии ??:

I I Vi2 h2 Vi2. (2.16) c1 h (V ) d c2 h i=1 i= Обозначив через i (x, y, z) кусочно-линейную функцию, равную 1 в i-ом узле сетки и нулю во всех остальных, произвольную кусочно-линейную функ цию можно записать в виде I h (2.17) V (x, y, z) = Vi i (x, y, z).

i= При подстановке этого выражения в (2.13) получаем I h W (V ) = Vi Vj (x, y, z)grad i (x, y, z) · grad j (x, y, z)d i,j= I 2 Vi i (x, y, z)q(x, y, z) + i= +(x, y, z)grad i(x, y, z) · grad (x, y, z) d. (2.18) Обозначим значения интергралов через Lij и qi. Тогда I I h (2.19) W (V ) = Lij Vi Vj 2 Vi qi.

i,j=1 i= Хотя интегрирование в (2.18) производится по всей области, фактически grad vi(x, y, z) = 0 только в тетраэдрах, одной из вершин которых является iтый узел сетки. Поэтому отличны от нуля только Lij для соседних узлов i, j.

Поскольку внутри каждого тетраэдра vi(x, y, z) является линейной функ цией, ее градиент есть постоянный вектор, который может быть вычислен через координаты вершин.

Пронумеруем вершины тетраэдра как показано на рис. 2.1.

j                 i Рис. 2.1. Нумерация вершин тетраэдра.

Обозначим через Ri, Rj, R1, R2 радиус-векторы вершин.

Найдем единичные нормали к граням, противоположным вершинам i, j:

ni = [(R2 Rj ) (R1 Rj )]/|[(R2 Rj ) (R1 Rj )]| (2.20) nj = [(R1 Ri ) (R2 Ri )]/|[(R1 Ri ) (R2 Ri )]|, и соответствующие высоты:

Hi = ((Ri Rj ), ni)ni (2.21) Hj = ((Rj Ri ), nj )nj.

Поскольку vi(x, y, z) = 0 на грани 1,2,j и vj (x, y, z) = 0 на грани i,1,2, в рассматриваемом тетраэдре grad vi(x, y, z) = Hi /Hi (2.22) grad vj (x, y, z) = Hj /Hj.

Поэтому скалярное произведение в интеграле (2.18) равно grad vi(x, y, z) · grad vi (x, y, z) = (Hi, Hj )/(Hi2Hj ).

(2.23) После вынесения этой константы из-под интеграла для тетраэдра остается (2.24) (x, y, z)d, то есть среднее значение, умноженное на объем тетраэдра (2.25) tetr tetr.

Объем тетраэдара можно выразить как (2.26) tetr = ((Rj Ri ), [(R1 Ri ) (R2 Rj )])) /6.

Таким образом, первый интеграл в (2.18) по тетраэдру равен:

tetr tetr (Hi, Hj )/(Hi2Hj ) (2.27) Просуммировав эти значения по всем тетраэдрам, прилегающим к отрезку ij, получаем значение Lij.

Параметры узлов i и j входят в выражение (2.27) симметрично, как и функции vi(x, y, z), vj (x, y, z) в вычисляемом интеграле, вклады рассматри ваемого тетраэдра в Lij и Lji совпадают. Этим обеспечивается Lij = Lji, то есть симметрия матрицы L.

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

Можно заметить, что если функция q(x, y, z) постоянна в тетраэдре, то qi = vi(x, y, z)q(x, y, z)d + grad vi (x, y, z) · grad (x, y, z)d = = qtetr tetr /4 + tetr Hi · grad (x, y, z)tetr/Hi2. (2.28) Коль скоро все Lij и qi есть известные числа, функционал энергии (2.19) есть квадратичная форма I числовых параметров Vi. Условием минимально сти его значения является равенство нулю всех производных по этим свобод ным параметрам, то есть система линейных алгебраических уравнений I (2.29) Lij Vj = qi, i = 1,..., I.

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

Каждое из линейных алгебраических уравнений (2.29) связывает значе ние Vi в данном узле и значения во всех соседних узлах. Как уже отмечалось, матрица L получается симметричной. Матрица получается положительно определенной, поскольку кусочно-линейные функции принадлежат энерге тическому пространству, в котором выполняется неравенство (2.15) и первое неравенство (2.16).

I I Vi2. (2.30) Vi Lij Vj c1 c1 h i,j=1 i= Оценки собственных чисел этой матрицы только множителями 1 и 2 из ограничений (2.5) отличаются от оценок в случае уравнения Пуассона.

2.5. Построение сетки.

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

При моделировании атмосферы в целом приходится решать задачу в об ласти, близкой к шаровому слою. Поэтому сетку строим из лучей, проведен ных из центра Земли через некоторые точки на сфере. Пример разбиения по высоте на 16 слоев показан на Рис. 2.2. На этом же рисунке показано вы сотное распределение атмосферной проводимости, из обзора [96], которое мы h, км 1014 1012 1010 108, См/м Рис. 2.2. Высотное распределение проводимости в атмосфере из [96]. Справа показано разбиение на 16 слоев по высоте.


используем в настоящей главе для демонстрации возможностей созданного численного метода.

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

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

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

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

Рис. 2.4. Разбиение шестигранника на 5 тетраэдров.

Простейшим является разбиение на 5 тетраэдров, показанное на рис. 2.4.

Четыре жирные точки являются вершинами тетраэдров, имеющих по три треугольных грани, лежащих на границе рассматриваемой шестигранной ячейки сетки. У пятого тетраэдра на границе лежат только ребра. Посколь ку в соседней ячейке должно быть согласованное разбиение на тетраэдры, в ней четырехугольники будут разрезаны на треугольники другими диагона лями. В результате в сетке получатся узлы двух видов. Если они находятся внутри области, то узлы, отмеченные жирными точками, имеют 6 соседних узлов, а остальные – 18. Поэтому уравнения (2.19), соответствующие этим уз лам, будут иметь разное количество слагаемых. Более существенным неудоб ством является связанное с выделением одной из диагоналей ограничение на объединение подобластей с такими разбиениями. Например, нельзя замкнуть такой криволинейный прямоугольник в кольцо при четном числе ячеек.

7 3 4 Рис. 2.5. Разбиение шестигранника на пирамиды по 4 тетраэдра в каждой.

Чтобы диагонали были равноправными мы выбираем некоторый центр шестигранника и строим 6 четырехугольных пирамид, основанием которых являются грани исходного шестигранника. В четырехугольном основания каждой пирамиды также выбираем некоторый центр и строим 4 треуголь ника, как это показано для одной грани ячейки сетки на рис. 2.5. Всего в ячейке получается 24 тетраэдра.

При объединении таких шестигранников в более сложные области един ственным ограничением является совпадения числа ячеек на общих гранях.

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

Если рассматривать значения V h (x, y, z) в построенных дополнительных узлах как независимые параметры, в системе уравнений (2.19) будет раз личное количество слагаемых в уравнениях, соответствующих узлам разного типа, что не удобно. Поэтому мы получаем эти значения с помощью интер поляции из восьми вершин шестигранной ячейки или из четырех вершин соответствующих грани. Обозначим значения Vi в этих основных K узлах как Uk. Все эти K узлов лежат внутри области, поскольку потенциал в гра ничных узлах равен нулю в силу однородных граничных условий (2.11, 2.11).

Запишем совокупность интерполяционных выражений, добавив к ним тож дества для K узлов, являющихся вершинами шестигранной ячейки сетки, в виде K (2.31) Vi = aik Uk, k= где при каждом i от нуля отличны восемь, четыре или одно значение aik.

В дальнейшем нам понадобится неравенство I K Vi2 (2.32) Uk, i=1 k= очевидное, поскольку среди Vi есть все Uk, и противоположная оценка, для которой существенен способ интерполяции (2.31). Пусть все коэффициенты интерполяции превосходят соответствующие средним арифметическим зна чения 1/4 и 1/8 не более, чем в 1 раз. Для положительности остальных трех или семи коэффициентов в интерполяционной формуле должно быть 4 или 8. Фактически мы используем 1. Для узлов типа 9 и типа 15 на рис. 2. (2.33) V9 = a1 U1 + a2 U2 + a3 U3 + a4 U4, V15 = a1 U1 +... + a8 U8, где для наглядности использована нумерация узлов на этом рисунке. Возве дя в квадрат, и используя предположенную ограниченность коэффициентов усреднения, получаем неравенства V92 2 (U1 + U2 + U3 + U4 )/4, 2 2 2 (2.34) V15 2(U1 +... + U8 )/8.

2 2 (2.35) Поскольку каждое значение Uk входит в оценки 12 значений типа 9 (2.34) и 8 значений типа 15 (2.35), после суммирования этих неравенств с добавлением неравенств вида V12 2 U1, следующих из тождеств вида V1 = U1, получаем дополнительное к (2.32) неравенство I K Vi2 2 (2.36) 5 Uk.

i=1 k= Функционал энергии (2.19) можно переписать в виде I K K I K h W (V ) = Lij aik Uk ajp Up 2 qi aik Uk = p= i,j=1 i= k=1 k= K I K I = aik Lij ajp Uk Up 2 qiaik Uk = i,j=1 i= k,p=1 k= K K (2.37) = LkpUk Up 2 qk U k.

k,p=1 k= Поскольку теперь у нас только K свободных параметров, минимизация функционала энергии соответствует K равенствам нулю производных по ним, то есть получается системе K уравнений для K значений V h во внутренних узлах:

K (2.38) LkpUp = qk, k = 1,..., K.

p= Эта система имеет тот же вид, что и (2.29), но меньшее число неизвестных и уравнений. Ее коэффициенты и правые части получаются суммированием Lij и qi с коэффициентами интерполяции aik по правилам (2.37).

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

Для всякого набора значений Up в узлах сетки с шестигранными ячейка ми по правилам интерполяции (2.31) получается набор значений Vi в сетке, состоящей из тетраэдров. В силу определения Lkp равны квадратичные фор мы I K (2.39) Lij Vi Vj = Lkp Uk Up, i,j=1 k,p= и по построению они равны квадратичной части интеграла энергии (2.18), которая также может быть записана как квадратичная форма (2.14) для кусочно-линейной функции V h. В силу (2.30) из (2.39) получаем:

K I Vi2. (2.40) LkpUk Up c1 c1 h i= k,p= С учетом (2.32) можно усилить это неравенство и получить положитель ную определенность матрицы L:

K K 3 (2.41) LkpUk Up c1c1 h Uk.

k,p=1 k= На первый взгляд, это неравенство с тем же коэффициентом, что и нера венство (2.30) для матрицы L, можно было бы получить проще. Интерполя ционные формулы (2.31) можно трактовать как переход из Iмерного про странства Vi в Kмерное подпространство, в котором независимыми явля ются K параметров Uk. Минимальное собственное число в подпространстве не может быть меньше, чем во всем пространстве, и также максимальное собственное число не может возрасти. Однако это справедливо, только если нормировка в подпространстве - та же, что и в исходном пространстве. Ис пользуемые нами нормы, квадраты которых равны суммам квадратов Vi и Uk, этого свойства, очевидно, не имеют. Они связаны лишь неравенствами (2.32), (2.36), и предельные значения могут реализоваться довольно точно.

Например, если все Uk равны 1, то при любом приемлемом правиле интер поляции все Vi равны 1 за исключением ячеек сетки около границ области, где получаются меньшие значения из-за усреднения с нулевыми значениями в граничных узлах. Поскольку для таких сеток (2.42) I 5K, для этого случая получаем I K Vi2 (2.43) 5 Uk.

i=1 k= Если же значения Uk = ±1, и знак меняется при переходе в соседний узел, то интерполированные значения Vi близки к нулю, кроме приграничных ячеек, а в K общих узлах значения Vi и Uk совпадают. Для этого случая получаем I K Vi2 Uk.

i=1 k= Из (2.30) действительно следует (2.41), а максимальное собственное число матрицы L может быть больше максимального собственного числа матрицы L не более чем в 5 раз. Это слишком грубая оценка, поскольку возрастание нормы в 5 раз (2.43) происходит для наборов значений Uk, соответствующих не большим, а малым собственным числам.

Просуммируем все уравнения (2.29), умножив их на Vi. Получаем I I (2.44) Lij Vi Vj = Vi qi.

i,j=1 i= Неравенство (2.30) позволяет оценить левую часть снизу, а стоящее в пра вой части скалярное произведение не превосходит произведения длин векто ров. Получаем неравенство I I I Vi2 Vi2 c1c1 h qi, i=1 i=1 i= из которого следует ограниченность решения системы уравнений (2.29) I I 1 qi Vi2 (2.45).

h c1c i=1 i= Для решения (2.38) из (2.41) аналогично получаем K K 1 qk 2 (2.46) Uk.

h c1 c k=1 k= Напомним, что qi и qk есть интегралы плотности источника тока q(x, y, z) по определенным окрестностям i-того и k-того узлов сетки. Поэтому для функции q(x, y, z), модуль которой ограничен, эти неравенства оказывают ся фактически независимыми от шага сетки.

Матрица L, как и L, является симметричной, поскольку правило (2.37) можно трактовать как умножение симметричной матрицы L на прямоуголь ную матрицу a справа и на сопряженную к ней слева.

Система линейных алгебраических уравнений (2.38) называется вариационно-разностной схемой. Ее коэффициенты и правые части вы числяем по приведенным выше правилам. В параграфе 2.7. приведем значения для простейшей сетки с одинаковыми кубическими ячейками.

Используемые нами сетки называются блочно-структурированными.

Их часто используют в сочетании с декомпозицией области, поскольку каждый блок соответствует своей подобласти, и это разбиение может быть рассмот рено как декомпозиция на непересекающиеся подобласти [15]. Декомпозиция особенно эффективна, когда есть специальные методы решения соответству ющих блоков системы сеточных уравнений (2.38), например, быстрое пре образование Фурье для каждой подобласти. При использовании достаточно произвольных криволинейных сеток эффективность декомпозиции обеспечи вается более сложными методами, включающими выбор граничных условий на внутренних границах, получающихся при разрезании, [15]. В некотором смысле можно сказать, что мы используем декомпозицию с условиями Дири хле. Однако можно трактовать наши действия просто как порядок нумерации неизвестных, описанный в приложении B B.6.. Для ускорения сходимости мы используем многосеточный метод.

2.6. Многосеточный метод.

Вариационная формулировка задачи позволяет обычным образом построить многосеточный метод. Мы следуем оригинальной версии метода, созданно го Р.П. Федоренко [44], с единственным отличием: вместо объединения трех ячеек по каждому направлению, мы стоим крупные сетки, исключая каж дый второй узел сетки в каждом направлении, пока не останется всего два шага сетки в одном из направлений, то есть на самой крупной сетке будет всего один внутренний узел на соответствующих сеточных линиях. Отметим, что для других краевых задач целесообразно провести еще одно укрупнение сетки, когда внутренних узлов не останется, а для задачи Дирихле на такой сетке уже не будет неизвестных. Использование для преобразования правых частей при укрупнении сетки оператора, сопряженного к оператору интерпо ляции, автоматически сохраняет симметрию матриц на крупных сетках.

Положительная определенность матриц на укрупненных сетках может быть доказана точно так же, как и (2.41), поскольку процедура интерполяции отличается от (2.31) только количеством членов и соотношением I 8K вме сто (2.42). Для такого важного параметра, как число обусловленности мат рицы, которое для симметричных матриц равно отношению максимального собственного числа к минимальному, доказывается невозможность возраста ния более, чем в 8 раз при укрупнении сетки вдвое. Эта оценка слишком груба, на самом деле, число обусловленности уменьшается примерно вчетве ро, что наряду с уменьшением числа неизвестных примерно в восемь раз при каждом укрупнении вдвое и обеспечивает эффективность многосеточного ме тода. Для классической разностной схемы на равномерной прямоугольной сетке для уравнения Пуассона это доказано в статье [44].

Многосеточный метод и способ его программной реализации описаны в Приложении A.2.3. Здесь отметим лишь, что на основной и вспомогатель ных сетках сглаженные решения строим с помощью 5 10 итераций метода верхней релаксации c параметром около 1.5.

В соответствии с теорией [44] для решения системы уравнений (2.38) тре буется порядка K lnK арифметических операций, где K – число узлов основ ной сетки. Иначе говоря, требуется порядка lnK многосеточных итераций.

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

2.7. Сетка с кубическими ячейками.

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

В тетраэдре с вершинами 1,9,4,15, показанном на рис. 2. V h V h V h V1 V4 V9 (V1 + V4 )/2 V15 V (2.47) =, =, =.

x h y h/2 z h/ Вклад этого тетраэдра в энергию (квадратичную часть функционала энергии (2.18)) равен h3 (V1 V4 )2 (V9 (V1 + V4 )/2)2 (V15 V9 ) (2.48) + +.

h2 h2 /4 h2 / Используя аналогичные выражения для всех 24-х тетраэдров, получаем вклад в энергию всей ячейки сетки:

h (Vi Vj )2, (2.49) 16 i,j= где для простоты записи включены и 8 тождественно нулевых членов при i = j. Во избежание недоразумений подчеркнем, что эти индексы i, j введены в одной ячейке сетки и не имеют отношения к нумерации узлов всей сетки.

Это выражение обращается в нуль, если V h (x, y, z) есть константа, и несложно показать, что для линейной функции V h (x, y, z) = V 0 Ex x Ey y Ez z, 0 0 (2.50) 0 0 где вектор (Ex, Ey, Ez ) есть напряженность электрического поля, однородно го в рассматриваемой ячейке сетки, энергия (2.49) равна h3 (Ex )2 + (Ey )2 + (Ez )2.

0 0 (2.51) Совпадение с точным значением говорит об отсутствии погрешности дис кретизации на линейных функциях.

Чтобы перейти к суммированию по ячейкам сетки, естественно использо вать нумерацию узлов сетки по трем направлениям. Пусть на рис. 2.5 узлы на левой грани имеют общий номер n, n + 1 – на правой, m – на ближней грани, m + 1 – на дальней, l – на нижней, l + 1 – на верхней. При этом нет смысла использовать обозначение U вместо V в этих узлах, как это было необходимо в параграфе 2.5. Интервалы нумеруем полуцелыми индексами.

Тогда изображенная ячейка сетки на рис. 2.5 имеет номер n + 1/2, m + 1/2, l + 1/2, и ее вклад в энергию (2.49) может быть записан в виде:

h (Vnml Vn+,m+,l+ )2. (2.52) n+1/2,m+1/2,l+1/,,= Вычисляя производную по Vnml, получаем вклад ячейки в уравнение (2.38) с номером k, соответствующим тройному номеру n, m, l:

h (2.53) n+1/2,m+1/2,l+1/2 (Vnml Vn+,m+,l+ ).

,,= Для пары узлов n,m,l и n + 1, m, l коэффициент Lkp складывается из четырех ячеек, прилегающих к соединяющих их ребру:

h L(n,m,l),(n+1,m,l) = (n+1/2,m+1/2,l+1/2 + n+1/2,m1/2,l+1/2 + (2.54) +n+1/2,m+1/2,l1/2 + n+1/2,m1/2,l1/2), или h, если проводимость одинакова в этих ячейках.

Для лежащих на одной грани узлов n, m, l и n + 1, m + 1, l суммируются вклады двух ячеек, прилегающих к грани l:

h (2.55) L(n,m,l),(n+1,m+1,l) = (n+1/2,m+1/2,l+1/2 + n+1/2,m+1/2,l1/2), или h/2 при const.

Пары узлов, являющиеся концами одной из главных диагоналей ячейки сетки, связаны только внутри этой ячейки, и h (2.56) L(n,m,l),(n+1,m+1,l+1) = n+1/2,m+1/2,l+1/2.

Аналогичные выражения получаются для всех пар узлов, и в силу (2.53) их суммированием и сменой знака суммы получаем коэффициент L(n,m,l),(n,m,l), соответствующий главной диагонали матрицы в записи (2.38).

Приведем уравнение (2.38), соответствующее узлу n, m, l, только при = const в окружающих этот узел ячейках сетки. Подразумевая, что узел не лежит на границе сеточной области, получаем:

h 14 Vn,m,l Vn+1,m,l Vn1,m,l Vn,m+1,l Vn,m1,l Vn,m,l+1 Vn,m,l 1 1 1 1 Vn,m+1,l1 Vn,m1,l+1 Vn,m+1,l+1 Vn,m1,l1 Vn+1,m1,l 2 2 2 2 1 1 1 1 Vn1,m+1,l Vn+1,m+1,l Vn1,m1,l Vn+1,m,l1 Vn1,m,l+ 2 2 2 2 1 1 1 Vn+1,m,l+1 Vn1,m,l1 Vn1,m1,l1 Vn1,m1,l+ 2 2 4 1 1 1 Vn1,m+1,l1 Vn1,m+1,l+1 Vn+1,m1,l1 Vn+1,m1,l+ 4 4 4 1 Vn+1,m+1,l1 Vn+1,m+1,l+1 = qnml. (2.57) 4 Как видим, в это уравнение вовлечены значения искомой кусочно линейной функции V h (x, y, z) во всех 27 узлах сетки, включая сам централь ный узел n, m, l.

Проанализируем полученную вариационно-разностную схему (2.38) с точ ки зрения теории разностных схем [2]. Предварительно напомним, что при простейшей аппроксимации производных получается разностная схема Vn+1,m,l + Vn1,m,l + Vn,m+1,l + Vn,m1,l + Vn,m,l+1 + Vn,m,l1 6Vn,m,l = qnml, h (2.58) имеющая второй порядок аппроксимации [2]. Ее очевидное преимущество по сравнению с нашей схемой (2.57) в меньшем количестве слагаемых, однако в случае неравномерных сеток оно исчезает. Это не доказывает отсутствие дру гих достоинств, которые могут быть связаны с лучшей аппроксимацией диф ференциального уравнения, но мы используем вариационно-разностный под ход, поскольку он представляется более простым для реализации на нерав номерных сетках.

Предполагаем, что узел n, m, l является началом декартовых координат, направления осей которых соответствуют увеличению индексов n, m, l. Под ставим в уравнение (2.57) ряды Тейлора узловых значений искомой функции V (x, y, z), которые для узлов трех типов аналогичны выражениям:

h2 2V h3 3V V + O(h4 ), (2.59) Vn+1,m,l = Vn,m,l + h + + 2 x2 6 x x h2 2V 2V 2V V V Vn+1,m+1,l = Vn,m,l + h + + +2 + + 2 x2 y x y xy (2.60) h3 3V 3V 3V 3V + +3 2 +3 + + O(h ), 6 x3 xy 2 y x y V V V Vn+1,m+1,l+1 = Vn,m,l + h + + + x y z h2 2V 2V 2V 2V 2V 2V + + + +2 +2 +2 + 2 x2 y 2 z 2 xy xz yz (2.61) h3 3V 3V 3V 3V 3V + + + +3 2 +3 + 6 x3 y 3 z 3 xy x y 3V 3V 3V 3V 3V + O(h4 ).

+3 2 + 3 2 + +3 2 + 3 + x z x z y z yz xyz Последние слагаемые в каждом из выражений являются остаточными членами в форме Лагранжа и вычисляются в некоторых точках внутри соот ветствующих интервалов (ребер рассматриваемой кубической ячейки сетки), тогда как производные до третьего порядка – в узле n, m, l. Например, в (2.59) остаточный член равен h4 4V O(h4 ) = (2.62), 24 x4 (x,0,0) где x есть некоторое число в интервале 0 x h. В остальных выраже ниях это более сложные комбинации четвертых производных. В предположе нии ограниченности всех четвертых производных искомой функции V (x, y, z) константой M, остаточные члены в приведенных рядах Тейлора по модулю не превосходят величин Mh4 /24, 2Mh4 /3, 32Mh4 /3, соответственно.

Заметим, что при нулевой правой части, qnml = 0, значение Vnml есть среднее от значений в 26 соседних узлах. Коэффициенты усреднения равны 1/14 для 6 значений, 1/28 для 12 значений и 1/56 для 8 значений. Все они положительны, и в сумме равны 1. Поэтому справедлив принцип максимума:

Vnml не превосходит максимального значения V h в этих узлах. Как следствие, для однородной задачи, когда все числа qnml = 0, решение имеет максимум и минимум на границе сеточной области, и все Vnml = 0, поскольку V h = 0 на границе. Из этого, в частности следует единственность решения неоднород ной системы линейных алгебраических уравнений (2.38) в рассматриваемом частном случае, когда уравнения принимают вид (2.57).

2.8. Разностная схема.



Pages:   || 2 | 3 | 4 |
 





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

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