Вычислительные технологии
Том 14, № 6, 2009
Математическое моделирование магнитосферных явлений и их влияния на атмосферу Земли*
В. В. ДЕНИСЕНКО, Н. В. ЕРКАЕВ Учреждение Российской академии наук Институт вычислительного моделирования СО РАН, Красноярск, Россия e-mail: [email protected], [email protected]
Описаны математические модели процессов генерации электрического поля в магнитосфере Земли. Изучены колебания токового слоя магнитосферного хвоста, возникающие при наличии градиентов компонент магнитного поля. Предложена модель электропроводности атмосферы с замыканием токов через ионосферу.
Ключевые слова: магнитосфера, атмосфера, электрическое поле, математическая модель.
В Институте вычислительного моделирования СО РАН с самого его основания большое внимание уделяется исследованиям в области математического моделирования физических процессов, происходящих в околоземном космическом пространстве. Построены модели таких важнейших процессов как обтекание солнечным ветром области геомагнитного поля [1-3] и генерация электрического поля за счет движения плазмы в хвосте магнитосферы [4]. Созданы модели магнитосферного магнитного поля [5] и ионосферных электрических полей и токов [6]. В основном эти модели предназначены для описания крупномасштабных процессов, для которых применимы уравнения магнитной гидродинамики. При рассмотрении отдельных магнитосферных явлений или объектов используются возможные упрощения: квазистатические приближения, выделение пограничных слоев, переход к двумерным моделям. В большинстве случаев приходится строить решения численными методами и иногда создавать новые численные алгоритмы, такие как многосеточный метод для решения двумерных задач переноса в гиро-тропных и движущихся средах, основанный на энергетических формулировках задач [7, 8].
Состояние высокоширотной ионосферы существенно зависит от физических процессов, протекающих в хвосте магнитосферы. Поэтому математическое моделирование нестационарных процессов в хвосте магнитосферы является важной и актуальной задачей. В работах [9-11] представлена новая модель изгибных магнитогидродинамических колебаний токового слоя магнитосферного хвоста, связанных с градиентами компонент магнитного поля. Такого типа колебания, получившие название флэппинг-колебаний, неоднократно регистрировались космическими аппаратами CLUSTER.
В рамках построенной нами модели [9] проведены детальные исследования модели флэппинг-волн, поведение которых характеризуется произведением градиентов тангенциальной и нормальной компонент (по отношению к плоскости токового слоя) магнит-
* Работа выполнена при поддержке РФФИ (гранты № 07-05-00135 и 09-06-91000) и Программы РАН № 16.3.
© ИВТ СО РАН, 2009.
ного поля магнитосферного хвоста. Флэппинг-волны распространяются вдоль токового слоя в направлении флангов перпендикулярно магнитному полю. Скорость распространения этих волн монотонно убывает с уменьшением длины волны. На рис. 1 показаны частота, групповая и фазовая скорости в зависимости от волнового числа. Частота колебаний монотонно растет с увеличением волнового числа и асимптотически стремится к величине ш>/, определяемой формулой
_ I 1 B * dBz
ш _ у 4пр X ~зХ'
где р — плотность, B* — максимальная тангенциальная компонента магнитного поля (по отношению к токовому слою), Bz — нормальная компонента невозмущенного магнитного поля, А — полутолщина токового слоя, ось x направлена вдоль слоя.
Результаты математического моделирования сравнивались с имеющимися экспериментальными данными спутников CLUSTER. Полученные теоретические значения характерных частот и скоростей распространения флэппинг-колебаний токового слоя магнитосферного хвоста хорошо согласуются с существующими эмпирическими оценками, основанными на данных наблюдений. Расчетные значения скорости распространения изгибных колебаний поперек плазменного слоя магнитосферного хвоста находят-
(1)
Рис. 1. Частота, групповая и фазовая скорость флэппинг-колебаний в зависимости от волнового числа для изгибной (kink — сплошная линия) и симметричной (sausage — штриховая линия) мод колебаний
ся в пределах нескольких десятков километров в секунду. Показано, что токовый слой магнитосферного хвоста становится неустойчивым относительно изгибных деформаций в области локального уменьшения толщины слоя. В этом случае возмущения токового слоя никуда не распространяются, а лишь экспоненциально растут по амплитуде, достигая нелинейной стадии. Такие локальные перетяжки в токовом слое обычно предваряют появление импульсов магнитного пересоединения в хвосте магнитосферы. Нестационарные изгибные деформации плазменного слоя создают альфвеновские волны, которые распространяются вдоль магнитных силовых линий и приносят электромагнитные возмущения в ионосферу Земли. Для описания этих волн разработана соответствующая модель магнитной струны [12].
Флэппинг-колебания токового слоя инициируются ускоренными потоками плазмы, возникающими вследствие импульсных процессов пересоединения магнитных полей в хвосте магнитосферы Земли. Существуют многочисленные экспериментальные подтверждения существования таких потоков. Скорость этих потоков может достигать нескольких сотен километров в секунду. Движущиеся к Земле с большой скоростью плазменные образования возбуждают магнитогидродинамические колебания плазменного слоя большой амплитуды и порождают флэппинг-колебания. В работе [13] представлена аналитико-численная модель пересоединения магнитных полей с учетом эффектов Холла. Показано, что эффекты Холла играют ключевую роль в процессах магнитного пересоединения. Проведено детальное сравнение результатов аналитического магнитогидродинамического моделирования с результатами двумерного численного эксперимента методом "частиц в ячейках", основанным на прямых расчетах движения частиц совместно с решением уравнений Максвелла. Сравнение показало хорошее согласие распределений токов и магнитных полей, полученных на основе обеих моделей.
Принято считать, что магнитное пересоединение инициируется кинетическими плазменными неустойчивостями в токовом слое. Развитие неустойчивостей существенно зависит от вида начальной функции распределения частиц по скоростям. Обычно в теориях неустойчивостей в качестве начальных функции рассматривались распределения Максвелла, являющиеся равновесными для столкновительной плазмы. Однако в последние годы активно развивается теория так называемых каппа-распределений частиц по скоростям, которые более адекватны космической плазме, чем распределения Максвелла. В связи с этим в работе [14] была рассмотрена двухпотоковая плазменная неустойчивость в предположении каппа-функции распределения протонов по скоростям. Изучение влияния плазменных неустойчивостей на процессы магнитного пересоединения, формирование ускоренных потоков плазмы и возбуждение флэппинг-колебаний является предметом дальнейших исследований.
В последние годы в соответствии с общемировой тенденцией мы уделяем большее внимание процессам взаимодействия ионосферы с атмосферой и литосферой. Прикладная направленность моделирования этих процессов обусловлена желанием использовать космические средства для обнаружения предвестников землетрясений. Имеются экспериментальные данные об изменениях электрического поля в приземной атмосфере накануне землетрясений. Покрыть Землю сетью наземных датчиков проблематично, поэтому возник вопрос о возможности судить об этих полях на основе спутниковых измерений в верхних слоях ионосферы.
В ионосфере среда гиротропна, поскольку вектор магнитного поля задает выделенное направление вращения, и проводимость является тензором с. Закон Ома может быть записан отдельно для направления, параллельного вектору магнитной индук-
ции B, и для перпендикулярных компонент векторов плотности тока j и напряженности электрического поля E:
j\\ = a\\E\\, (2)
j± = apE± - an [E± x B]/B, (3)
где ap, an, a\\ — педерсеновская, холловская и продольная проводимости. Для вычисления значений этих компонент тензора проводимости создана модель, представленная в работе [15]. В качестве ее элементов использованы известные, свободно распространяемые эмпирические модели, которые дают пространственные распределения электронной и ионных концентраций IRI, газового состава и температуры MSISE и геомагнитного поля IGRF-10. Реализующие их программы на языке Фортран были взяты на сайте НАСА [16]. На рис. 2 показаны полученные высотные распределения компонент тензора проводимости, типичные для среднеширотной ионосферы выше 90 км.
В атмосфере среда изотропна, поэтому проводимость является скаляром: ap = a\\, an = 0. Ее высотный профиль построен на рис. 2 ниже 60 км в соответствии с эмпирической моделью [17]. В области высот от 60 до 90 км использовались интерполяционные формулы, позволяющие гладко сшить атмосферную и ионосферную модели.
Как видно из рис. 2, в ионосфере проводимость в направлении магнитного поля a\\ на несколько порядков превышает проводимости в перпендикулярных направлениях. Поэтому, как и во многих ионосферных моделях крупномасштабных электрических полей, мы приближенно считаем магнитные силовые линии эквипотенциальными. В таком
10-15 10-1О 10-5 а> См/м 10_1Б 10-Ю 10-5 См/м
Рис. 2. Распределения компонент тензора; проводимости в атмосфере и в ионосфере
случае получается двумерная модель, в которой каждая точка представляет целую силовую линию. Такая ионосфера эквивалентна тонкой проводящей пленке с двумерным законом Ома
л=(£ Iй) Е (4)
где интегральные педерсеновская и холловская проводимости получаются из локальных интегрированием вдоль силовой линии. Если магнитное поле не вертикально, закон Ома сохраняет этот вид в специальных координатах [15].
Используемая упрощенная модель ионосферы позволяет сформулировать для верхней границы атмосферы условие, соответствующее закону сохранения заряда: приходящий из атмосферы ток растекается по ионосфере. Если такой границей считать сферу радиуса т\, то
^ J = Зг\Г=Г1 . (5)
В отличие от этого естественного условия в известной модели, обзор развития которой дан в [18], без обоснований полагается
Зг \Г=Г1 = 0, (6)
и Т1 соответствует высоте 90 км над землей. Точнее говоря, поставлено эквивалентное условие
д Ф
дг
= 0,
Г=Г 1
где Ф — электрический потенциал такой, что напряженность электрического поля
Е = -УФ. (7)
Сравнение с (5) показывает, что условие (6) соответствует нулевым токам в ионосфере, т. е. ионосфере — изолятору, что явно противоречит обратному соотношению проводимостей ионосферы и атмосферы, показанному на рис. 2.
В обсуждаемой модели было предсказано проникновение в ионосферу электрических полей напряженностью до 1 мВ/м, которые сравнимы с обычно существующими полями иного происхождения.
В другой известной модели [19] предложено условие
Фи, =0
и в качестве граничной сферы радиуса т2 выбрана поверхность, лежащая внутри ионосферы на высоте более 150 км над землей. Это условие в силу (7) эквивалентно обращению в нуль касательных к границе компонент электрического поля и описывает границу с идеальным проводником, т. е. проводимость ионосферы выше 150 км полагается бесконечной. Авторами [19] адекватность такого приближения не проанализирована. Наши оценки показывают, что для замыкания полученного ими тока из атмосферы в ионосферу при полученных в области выше 100 км горизонтальных электрических полях проводимость ночной ионосферы должна быть в сотни раз больше, чем она есть.
Таким образом, две известные модели проникновения электрического поля от поверхности Земли до ионосферы приходится считать неудовлетворительными.
В созданной нами модели [20] учтена проводимость ионосферы. При этом использовалась та же модель электропроводности атмосферы, т. е. закон сохранения заряда, который с учетом (2), (3), (7) в изотропной среде имеет вид
—^у (а grad Ф) = 0.
(8)
На нижней границе атмосферы во всех трех моделях задается вертикальная компонента электрического поля, для чего используются некоторые представления, основанные на результатах измерений:
дФ
дг
(9)
г=г 0
В локальных задачах сферичность не существенна, и задача формулируется для плоского слоя с дополнительным условием убывания решения на больших расстояниях.
При тех же наземных полях напряженностью до 100 В/м, что и в модели [18], мы получили на порядки меньшие поля — не более 10 мкВ/м в ночной ионосфере, и еще в сто раз меньшие — в дневной. Величина проникающего поля оказывается примерно обратно пропорциональной интегральной проводимости ионосферы, а днем она на два порядка повышается от ночного уровня.
В настоящее время модель [20] усовершенствована для учета реального распределения проводимости в атмосфере, приведенного в [17], тогда как ранее мы задавали экспоненциальный рост с высотой. Рассмотрена и более реалистичная модель распределения По сути эта функция описывает подземный генератор тока, поскольку нормальная компонента плотности тока аЕ<0(9,^) сохраняет свое значение ниже поверхности Земли. Для локальных явлений, естественно, оба полюса такого генератора должны быть в области, где эти явления происходят, и ни один из них не может быть вынесен на бесконечность. В такой усовершенствованной модели для ионосферы получены еще в несколько раз меньшие электрические поля, проникающие в ионосферу.
На рис. 3 приведены линии тока в атмосфере. Часть тока, замыкающегося через ионосферу, существенно убывает с уменьшением горизонтального размера участка земной поверхности, на котором возникает электрическое поле. В приведенном на рис. 2 случае до ионосферы доходит около 30 % тока, выходящего из-под земли. При увеличении горизонтального размера источника тока вчетверо уже почти весь ток (97 %) течет через атмосферу вверх, где замыкается через хороший проводник, которым является ионосфера. В последнем случае мы получили в ночной ионосфере ионосферные поля не более 2 мкВ/м.
80
г, км
-100
X, км
100
Рис. 3. Линии тока в атмосфере
Выделение полей полученного нами масштаба на фоне всегда существующих в ионосфере полей более 1 мВ/м невозможно. Таким образом, к сожалению, исключена возможность использования спутниковых измерений для оценки наземных электрических полей. Разумеется, это не относится к другим известным механизмам связи приземной атмосферы с ионосферой, которые можно было бы использовать для обнаружения предвестников землетрясений [17].
В рамках рассматриваемого прикладного направления необходим также анализ обратного явления — проникновения ионосферных электрических полей до поверхности Земли, поскольку эти поля могут объяснять часть измеряемых наземных вариаций и после их исключения упростится выделение предвестников землетрясений. Поэтому мы создали комплекс программ для численного решения трехмерной стационарной задачи электропроводности атмосферы [21]. Пространственное распределение электрического потенциала получается как решение смешанной краевой задачи для эллиптического уравнения. Основная сложность задачи обусловлена изменением проводимости в расчетной области в миллионы раз. Для решения этой задачи была построена вариационно-разностная схема, которая получается стандартным образом из условия минимальности значения функционала энергии по узловым значениям кусочно-линейных функций, используемых для аппроксимации энергетического пространства. В созданном алгоритме применяются блочно-структурированные сетки с такой же нумерацией узлов в каждом блоке, как и в параллелепипеде. Каждая ячейка сетки, имеющая восемь вершин, после определения центра ячейки и центров ее шести граней делится естественным способом на 24 тетраэдра. Для решения возникающей системы линейных алгебраических уравнений с симметричной положительно определенной матрицей применяем многосеточный метод Федоренко. Потенциал на верхней границе атмосферы задается как потенциал в ионосфере. Для этого мы использовали модели ионосферных электрических полей, ранее созданные для различных геомагнитных условий [6]. Земля приближенно рассматривалась как идеальный проводник, поскольку проводимости морской воды и влажной почвы существенно превышают проводимость воздуха в приземном слое. Комплекс программ позволяет также решать задачи с иными граничными условиями. Для задания пространственного распределения проводимости были использованы существующие эмпирические модели [17].
В результате моделирования показано, что ионосферные разности потенциалов до 100 кВ, характерные для магнитосферных суббурь, приводят к наземным вариациям вертикальной компоненты электрического поля до 60 В/м. Это весьма существенная поправка к наблюдаемому при хорошей погоде среднему полю около 130 В/м, которое соответствует постоянной разности потенциалов между землей и ионосферой, обусловленной глобальной грозовой активностью. Методическим результатом моделирования является обоснование обычного способа построения вертикального приземного электрического поля в рамках одномерной модели электропроводности атмосферы для полей с достаточно большим горизонтальным масштабом поля. Одномерное приближение можно использовать до высоты 20 км, если горизонтальный масштаб превышает 50 км, и до высоты 50 км при горизонтальном масштабе более 500 км. Масштаб мы понимаем в том же смысле, в котором для функции вт(ж) характерным является изменение аргумента на единицу. Для построения полей в верхней атмосфере решение существенно трехмерной задачи является необходимым, и созданный комплекс программ позволяет это делать. Кроме того, описанный комплекс дает возможность учесть рельеф местности и сложные распределения проводимости в пространстве.
Главным недостатком созданной модели является упрощенный учет слоя, обычно занимающего высоту 80-100 км, где проводимость уже не является изотропной, как в атмосфере, но еще нет подавляющего преобладания проводимости вдоль магнитного поля, характерного для расположенной выше данного слоя основной части ионосферы. Для этого слоя необходимо решать следующие уравнения электропроводности с несимметричным тензором проводимости (2), (3):
¿IV j = <<,
го1 Е = О, Зп |г = д.
(10)
Ненулевая правая часть < в законе сохранения заряда возникает, если учесть заданную скорость движения проводника V, поскольку в законе Ома (3) добавится член, соответствующий переходу в движущуюся систему отсчета проводника и < = — (о^ х В]). Для ионосферы этот эффект важен, так как именно ветрами в рамках динамотеории объясняются основные токи, наблюдаемые в спокойных геомагнитных условиях.
Ненулевая правая часть в уравнении индукции О возникает при решении квазистационарных задач, когда магнитное поле заданным образом изменяется со временем. В этом случае О = —дВ/д£.
Граничное условие (3) эквивалентно (9), поскольку отличается лишь умножением на строго положительную проводимость, но является более содержательным, так как именно нормальная компонента плотности тока j непрерывна при переходе за границу. Вторым основным условием из возможных краевых условий является задание касательных компонент напряженности электрического поля вместо (4):
Ег |г
gт
поскольку они тоже непрерывны при переходе за границу.
Чтобы избежать решения задач с несамосопряженными операторами, которые традиционно получаются для электрического потенциала, задачу (2), (3), (11) целесообразно решать в рамках энергетической формулировки [7]. В соответствии с энергетическим методом [22] в [7] был построен функционал энергии
ж р)=2
F р
F р
— + Р*О) dQ + Fg/aodГ,
(11)
который рассматривается на множестве пар гладких функций F, Р, удовлетворяющих условиям
Fdtt = 0, Рг|г = 0.
(12)
Квадратными скобками в (11) и ниже обозначена симметричная билинейная форма, определяющая энергетическое скалярное произведение:
grad и
го1 V
—2 о Бо*
о,
1
--Бо*
оо
F Р
--оБ
оо
Б
grad F го1 Р
\
+ V Р
dQ.
/
Здесь положительная константа a0 и симметричная равномерно положительно определенная матричная функция S используются для улучшения обусловленности возникающих при численном решении задачи систем линейных алгебраических уравнений.
Полученные в результате минимизации функционала энергии функции F, P позволяют построить решение исходной задачи по формуле
E = -—Sa*grad F + S rot P. Условием минимума является выполнение уравнений
div (--2aSa*gradF +--SrotP ) = Q/a0,
V aö a0 J
(13)
rot —Sa*grad F + S rot P^ = G.
Еще одно следствие минимальности значения функционала энергии —
div P = 0.
В пространстве функций, удовлетворяющих этому условию и (12), оператор системы (14) является симметричным, положительно определенным и сопряженно факторизо-ванным, что существенно упрощает численное решение задачи по сравнению с традиционной формулировкой задач с несамосопряженными операторами для электрического потенциала.
Для задания проводимости в рассматриваемом слое нижней ионосферы, где существенна гиротропия, может быть использована модель [15], основанная на известных эмпирических моделях пространственно-временных распределений параметров ионосферы. Модель электропроводности этого слоя нижней части ионосферы будет частью общей модели, поскольку целесообразно провести декомпозицию всей области на три части, в двух из которых возможны существенные упрощения: за счет изотропии — в атмосфере, и за счет подавляющего преобладания проводимости вдоль магнитного поля — в основной части ионосферы. Достаточно эффективные численные реализации этих двух упрощенных моделей нами уже созданы.
Атмосферная часть задачи существенно проще, поскольку тензор проводимости a становится изотропным. Достаточно его симметрии, чтобы задача расщепилась на две независимые задачи для F и P. При нулевой правой части G решением второй задачи является P = 0. Функция F становится электрическим потенциалом, и получающееся из (14) уравнение для F совпадает с (8). Соответственно и используемый для решения этой задачи энергетический метод, кратко описанный выше, является частным и существенно более простым случаем минимизации функционала энергии (11).
Задача для основной части ионосферы сводится к двумерной и для нее энергетический подход реализуется в виде многосеточного метода [8].
Список литературы
[1] Денисенко В.В., Еркаев Н.В., Китаев А.В., МАТВЕЕНКОВ И.Т. Математическое моделирование магнитосферных процессов. Новосибирск: Наука, 1992.
[2] Еркаев Н.В. Обтекание солнечным ветром магнитосферы. М.: Междуведомств. геофиз. комитет АН СССР, 1989. 127 с.
[3] Erkaey N.V., Mezentsev A.V., Biernat H.K. Influence of the interplanetary magnetic field on the solar wind flow about planetary obstacles // Space Science Rev. 2006. Vol. 122. P. 209-219.
[4] Денисенко В.В., Замай С.С., Китаев А.В. Влияние вязкого трения между солнечным ветром и плазменным слоем на генерацию электрического поля в магнитосфере // Геомагнетизм и аэрономия. 2003. Т. 43, № 6. С. 730-736.
[5] Denisenko V.V., Biernat H.K., Erkaey N.V., Semenov V.S. Mathematical model of magnetic field perturbations by currents in the Earth's magnetosphere // Planetary Radio Emissions VI / Eds. H.O. Rucker, W.S. Kurth, G. Mann. Vienna: Austrian Acad. Sci. Press, 2006. P. 309-316.
[6] Denisenko V.V., Zamay S.S. Electric field in the equatorial ionosphere // Planetary and Space Science. 1992. Vol. 40, N 7. P. 941-952.
[7] Денисенко В.В. Энергетический метод для трехмерных эллиптических уравнений с несимметричными тензорными коэффициентами // Сиб. мат. журн. 1997. Т. 38, № 6. С. 1267-1281.
[8] Денисенко В.В. Энергетические методы для эллиптических уравнений с несимметричными коэффициентами. Новосибирск: Изд-во СО РАН, 1995. 204 с.
[9] Erkaey N.V., Semenov V.S., Biernat H.K. Magnetic double-gradient instability and flapping waves in a current sheet // Phys. Rev. Lett. 2007. Vol. 99. P. 235003.
[10] Erkaey N.V., Semenov V.S., Biernat H.K. Magnetic double gradient mechanism for flapping oscillations of a current sheet // Geophys. Res. Lett. 2008. Vol. 35, N L02111. doi:10.1029/2007GL032277.
[11] Erkaey N.V., Semenov V.S., Kübyshkin I.V. et al. MHD aspect of current sheet oscillations related to magnetic field gradients // Ann. Geophys. 2009. Vol. 27. P. 417-425.
[12] Еркаев Н.В., Шайдуров А.В. Модель магнитной струны для расчета колебаний тонких магнитных трубок // Вычисл. технологии. 2006. Т. 11, № 3. С. 73-89.
[13] Semenov V., Korovinskiy D., Divin A. et al. Collisionless magnetic reconnection: analytical model and PIC simulation comparison // Ann. Geophys. 2009. Vol. 27. P. 905-911.
[14] Langmayr D., Biernat H.K., Erkaey N.V. Influence of kappa-distributed ions on the two-stream instability // Phys. Plasmas. 2005. Vol. 12. P. 102103.
[15] Denisenko V.V., Biernat H.K., Mezentsev A.V. et al. Modification of conductivity due to acceleration of the ionospheric medium // Ann. Geophys. 2008. Vol. 26. P. 2111-2130.
[16] Models Distribution and Staging Directory. National Space Science Data Center. NASA. http: nssdcftp.gsfc.nasa.gov
[17] Molchanov O., Hayakawa M. Seismo-electromagnetics and Related Phenomena: History and Latest Results. Tokyo: TERRAPUB, 2008.
[18] Pulinets S.A., Legen'ka A.D., Gaivoronskaya T.V., Depuev V.Kh. Main phenomenological features of ionospheric precursors of strong Earthquakes //J. Atmosph. and Solar-Terrestrial Phys. 2003. Vol. 65. P. 1337-1347.
[19] Grimalsky V.V., Hayakawa M., Ivchenko V.N. et al. Penetration of an electrostatic field from the lithosphere into the ionosphere and its effect on the D-region before earthquakes // Ibid. 2003. Vol. 65. P. 391-407.
[20] Denisenko V.V., Boudjada M.Y., Horn M. et al. Ionospheric conductivity effects on electrostatic field penetration into the ionosphere. // Natural Hazards and Earth System Sci. J. 2008. Vol. 8. P. 1009-1017.
[21] Денисенко В.В., Бычков В.В., Помозов Е.В. Расчет атмосферных электрических полей, проникающих из ионосферы // Солнечно-земная физика: Сб. науч. тр. / РАН. Сиб. отд-ние. Ин-т солнечно-земной физики. 2008. Т. 2, вып. 12. С. 281-283.
[22] Михлин С.Г. Вариационные методы в математической физике. М.: Гостехиздат. 1957.
Поступила в редакцию 9 ноября 2009 г.