Конечноэлементное решение стационарной системы уравнений Максвелла с разрывными коэффициентами тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

КРЕМЕР Игорь Альбертович

КОНЕЧНОЭЛЕМЕНТНОЕ РЕШЕНИЕ СТАЦИОНАРНОЙ СИСТЕМЫ УРАВНЕНИЙ МАКСВЕЛЛА С РАЗРЫВНЫМИ КОЭФФИЦИЕНТАМИ

01.01.07 - вычислительная математика

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

- 2 СЕН 2010

Новосибирск - 2010

004607642

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

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

доктор физико-математических наук, старший научный сотрудник Урев Михаил Вадимович

Официальные оппоненты: доктор физико-математических наук,

профессор Кабанихин Сергей Игоревич доктор физико-математических наук, профессор Федорук Михаил Петрович

Ведущая организация:

Учреждение Российской академии наук Институт нефтегазовой геологии и геофизики им. A.A. Трофимука Сибирского отделения

РАН

Защита состоится «27» октября 2010 года в 17 часов на заседании диссертационного совета Д 003.061.01 при Учреждении Российской академии наук Институте вычислительной математики и математической геофизики Сибирского отделения РАН, по адресу: 630090, г. Новосибирск, проспект академика Лаврентьева, д. 6.

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

Автореферат разослан «5» июля 2010 года.

Ученый секретарь

диссертационного совета к. ф.-м. н.

Рогазинский С. В.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

rot a rot X + cßX - F, divßX = 0.

Здесь параметры a,ß зависят от электромагнитных свойств среды. В качестве неизвестной функции X могут выступать напряженности электрического или магнитного поля, а также векторный магнитный потенциал. Константа с > 0 зависит от гармонической частоты или от величины временного шага.

Традиционные схемы численного решения задачи (1) строятся на использовании регуляризирующего свойства слагаемого cßX. Можно решать только первое уравнение системы (1), а второе уравнение является простым следствием соленоидальности правой части div F = 0. Поскольку правая часть F зачастую определяется численными методами, то в реальных ситуациях дивергентные ограничения решения могут нарушаться в результате накопления погрешностей вычислений. Для устанавливающихся полей на поздних временах такие эффекты разрушают численное решение [3, 4]. Другой аспект проблемы связан с наличием нетривиального ядра у оператора rot, куда входят градиенты скалярных функций. При стремлении константы с к нулю, левая часть первого уравнения системы (1) становится вырожденной. В связи с вышесказанным задачу (1) будем рассматривать в предельном, стационарном варианте с = 0.

Для преодоления описанных проблем применяются различные методы регуляризации. В работах М. Costabel, М. Dauge (2002), J. Bramble, Т. Kolev, J. Pasciak (2004) развивается "grad - div" регуляризация, когда к оператору rot rot добавляется слагаемое —V(div ßX). Близким по идеологии является разрывный метод Галеркина (DG). В работах I. Perugia, D. Schotzau, P. Monk (2001), F. Kikuchi (1987) DG применяется для решения уравнений Максвелла

Г

\

\

в частотной области. Наиболее естественным функциональным пространством для искомых решений является H(rot; Элементы данного векторного пространства обладают свойством непрерывности касательных компонент и допускают разрыв нормальных компонент. В работах J.C. Nedelec (1980, 1986) предлагаются векторные "edge" элементы для аппроксимации H(rot; В методах "grad — div" регуляризации данное пространство не используется, поскольку на функциях X £ H(rot; Q) дивергенция div ßX е H-1(Q) является распределением, и возникают сложности интерпретации слагаемого —V(divßX). Авторами G. Haase, М. Kuhn, U. Langer (2001) при некоторых ограничениях получена регуляризированная дискретная постановка для стационарной магнитной задачи основанная на "edge" элементах. Аналогичные конструкции использовались С. Greif, D. Schotzau (2007), Q. Hu, J. Zou (2004) для численного решения задач с седловой точкой с вырожденными диагональными блоками. За рамками данных работ остались вопросы обоснования полученных регуляризированных постановок и сходимость численных решений к точным. В работах V. Girault, Р.-А. Raviart (1986) исследуется способ регуляризации абстрактных задач с седловой точкой.

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

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

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

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

Прикладная ценность работы. Результаты исследований реализованы в виде численных модулей и включены в программу «MODEM 3D» [5], предназначенную для интерпретации данных 3D нестационарных зондирований.

Апробация результатов диссертации. Результаты работы докладывались на объединенном семинаре кафедры вычислительной математики НГУ и ИВМ и МГ СО РАН (г. Новосибирск, рук. д.ф.-м.н. профессор Ильин В.П., 2009, 2010), в Baker Atlas (г. Новосибирск, рук. д.т.н. Табаровский Л.А., 2007). Также на Российских и международных семинарах и конференциях. 6-й Международный геофизический научно-практический семинар Применение современных электроразведочных технологий при поисках месторождений полезных ископаемых (г. Санкт - Петербург, 2008). 5-ая Международная конференция и выставка «Недра - 2008» (г. Москва, 2008). Научно -практическая конференция Комплексирование геолого - геофизических методов при обосновании нефтегазопоисковых объектов на Сибирской платформе (г. Новосибирск, 2008). Международная конференция «Математические Методы в Геофизике - 2008» (г. Новосибирск, 2008). 2-ая Всероссийская конференция «Актуальные проблемы строительной отрасли» (г. Новосибирск, 2009). Международная конференция по вычислительной математике 1ССМ - 2009 (г. Новосибирск, 2009). 1-ая международная конференция «Актуальные проблемы электромагнитных зондирующих систем» (г. Киев, 2009).

Публикации по теме диссертации. Результаты диссертации опубликованы в 7 работах, в том числе 5 работ в журналах из перечня ВАК [1, 2, 5, 6, 7] и две работы [3, 4] - тезисы международных конференций. Список работ приведен в конце автореферата.

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

Личный вклад автора. Основные результаты работы получены автором. Достаточные условия теоремы 2.3. установлены совместно с М.В. Уревым. Постановку практической задачи, приведенной в параграфе 3 четвертой главы, сделал М.И. Эпов. Программная реализация «MODEM 3D» выполнена совместно с М.И. Ивановым. Представленные и выносимые на защиту результаты, полученные в совместных исследованиях, согласованы с соавторами.

КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

Во введении (первая глава) рассмотрены различные постановки квазистационарной системы уравнений Максвелла в проводящих неоднородных областях. Обосновывается необходимость изучения стационарной системы уравнений

1

rot -rotA = F,

fi UJ

- divaA = 0.

Здесь A - векторный магнитный потенциал, F - источник электромагнитного поля, а - удельная электрическая проводимость, ц - магнитная проницаемость. Объектом наших исследований является стационарная система уравнений (2). Предметом исследований данной работы является решение таких систем методом векторных конечных элементов.

Введение также содержит обзор литературы. Кратко формулируются полученные результаты. Описывается структура диссертации.

Вторая глава содержит три параграфа, в которых проводится исследование обобщенной регуляризированной стационарной задачи в неоднородных проводящих областях. В §1 приводится схема решения абстрактной задачи с седловой точкой. Пусть дана пара гильбертовых пространств V, Q и определены непрерывные билинейные формы а(\•): V х V -* К и b(v): V X Q -* JR. Пусть F 6 V', рассматривается следующая задача. Требуется найти такую пару функций (и, р) £ V х Q, что выполнено:

a(u, v) + b(y, р) = (F, v>, Vv е V, b(u,q) = 0, VqtQ-

Если обозначить V0 = {и 6 V: b(v,q) = 0, Vq 6 (?}, то система уравнений (3) распадается на две задачи. Первая задача ставится в подпространстве У0: найти такую функцию и Е V0, что выполнено:

a(u,v) = (F,v), VveV0. (4)

Вторая задача, очевидно, следует из (3) и (4). Найти р G Q, что:

¿0, р) = (F, v) - а{и, v), Vv £ V. (5)

Достаточными условиями разрешимости задач (4) и (5) являются V0 - эллиптичность билинейной формы а(у) и выполнение "inf — sup" условия для билинейной формы Ь(-,-):

3аа > 0, a(v, v) > аа\\v\\$, Vv е V0,

b (v, q)

3ab > 0, supper -.. ,. >ab\\q\\Q, VqtQ-

В рамках такой схемы решения системы (3) получается, что необходимо решать задачу (4) в подпространстве функций С точки зрения численного решения этой задачи, как минимум, возникает необходимость строить конечномерные аппроксимации подпространства У0. Эти построения могут оказаться достаточно сложной отдельной задачей. Для преодоления данного затруднения используется прием регуляризации. Во второе уравнение системы (3) вносится малое возмущение и формулируется приближенная задача. Найти такую пару (и£,р£) 6 V х Q, что выполнено:

а(и£, v) + b(v, р£) = (F, v), Vv е V, b{u£,q)-Ec(p£,q) = Q, Vq G Q. ()

Здесь s > 0 - малый параметр, c(v): Q х Q -* Е Q - эллиптичная непрерывная билинейная форма. Если определить операторы В 6 £(V, Q') и се £((2,(20:

(Bv,q) = b(v,q), (Cp,q) = c(p,q), VvEV, Vp.qeQ, то из второго уравнения системы (6) можно выразить р£ как:

р£=-С~1Ви£. (7)

е

Для неизвестной функции и£ получаем задачу. Найти и£ 6 V, что: 1

а(и£, v) +-{Bv, С~1Ви£) = (F,v), Vv € V. (8)

В отличие от (4), решение задачи (8) определяется во всем пространстве V. Для достаточно малых значений £ справедлива оценка:

\\uE-u\\v + \\p£-p\\Q<K£\\F\\r. (9)

В §2 стационарная система уравнений для векторного потенциала (2) формулируется как обобщенная задача седловой точкой. Пусть ПсМ3 открытый выпуклый многогранник с границей 3fi. Предполагаем, что удельная проводимость а и магнитная проницаемость /г являются постоянными в подобластях положительными скалярными функциями:

О < сг0 < а < отах, 0<цо<ц< цтах.

В качестве неизвестных функций в задаче (3) выступают векторный потенциал и 6 V и множитель Лагранжа р 6 Q. Здесь решение определяется в функциональных пространствах

V = H0(rof,a), <2=я01(п).

Билинейные формы определяются следующим образом:

a(u,v) = f -rotu-rotvdil, Vu.veV,

J a^

b(v,q) = - au-Vqdn, Vv 6 V, Vq e Q.

В этом случае V0 - является множеством функций удовлетворяющих калибровочному условию div а А = 0 в слабом смысле.

В следующих двух леммах и теореме для введенных билинейных форм проверяются достаточные условия разрешимости задач (4) и (5) а, следовательно, и исходной задачи (3). Используя специфику магнитостатической задачи, устанавливаются априорные оценки решения. Лемма 2.1. Форма К X V К. V0 - эллиптична:

3аа>0, a(v,v) > aJMIrot.n Vv е V0. Лемма 2.2. Форма b(v): V х Q -> R удовлетворяет "inf — sup" условию:

b(v,q)

Баь > 0, sup„ey тг-г-> ablklli,n Vq € Q.

Теорема 2.1. Решение задачи (3) (u,p) Е V х Q существует и единственно. На решении выполнены оценки устойчивости

IMUnS —НЯ1г. Ilplli.il <-НЯ1г-

«а аЪ

В §3 для сформулированной в предыдущем параграфе системы уравнений применяется прием регуляризации. Регуляризация производится с использованием билинейной формы с(у): (? х <3 -* К.

¿п

При выводе регуляризированной задачи учитывается тот факт, что удается выделить задачу для множителя Лагранжа - р. Поскольку с V то в первом уравнении системы (3) в качестве пробных функций можно использовать V = Уд. Верно тождество р) — ~с(р, ц), В этом случае требуется найти такую функцию р Е ф, что выполнено

с(р,я) = (<1юР,ч), Уде (2. (10)

Из второго уравнения системы (3) и (10) следует точное равенство для любого значения параметра регуляризации р > 0:

Ь(и,ч)-рс(р,ч) = -р{(ИуР,я), V<г (И)

Тождество (11) дает выражение, аналогичное (7) (искомые величины помечаем индексом /?):

1

рр = — С'1 Вир + Р. (12)

Первое уравнение системы (3) и (12) дает регуляризированную задачу для векторного потенциала. Найти функцию ир е У такую, что выполнено:

а(ир,у) + С'1 Вир) = (Р,у) - (Вг, С_1<21Р Уу б V. (13)

Теорема 2.2. Для любого р > 0, задача (13) имеет единственное решение ир е V.

Для исследования разрешимости задачи (13) используется вариант разложения Гельмгольца.

Лемма 2.3. Для каждой функции и е V имеет место разложение

у = ы + (14)

Здесь IV еУ0, д 6 <?. Кроме того, справедлива оценка:

М2г„(,п £ Ф™ п\\20Д + Н^НЗл). ' (15)

Основной результат второй главы формулируется в виде теоремы. Теорема 2.3. Задачи (13) и (12) эквивалентны- исходной задаче (3) для любого положительного значения параметра регуляризации р.

Таким образом, решение регуляризированной задачи (13) удовлетворяет калибровочному условию вне зависимости от соленоидальности правой части. Устанавливается, что билинейная форма ар{-,-)\ V х V К, /? > 0, соответствующая регуляризированной задаче

1

ар (и, V) = а(и, v) + - (Ву, С'1Ви), Vu,veV,

симметрична: а,р(и, V) = ар (у, и).

В третьей главе рассматривается дискретная регуляризированная магни-тостатическая задача для векторного потенциала. В §1 вводится семейство согласованных, регулярных и квазиравномерных разбиений области П: Т/г = I = 1, - ,N1^. Здесь /1 = тахкет/1кк, кк - диаметр тетраэдра К. Определяются конечномерные аппроксимирующие пространства с V и ?1,с Пространство Уп строится с использованием векторных элементов Неделека первого порядка, второго типа. Степени свободы таких элементов связаны со значениями касательных компонент векторных функций на ребрах тетраэдров. Векторные функции уп £ характеризуются тем, что на гранях тетраэдров имеют непрерывные касательные компоненты, в то время как нормальные компоненты могут быть разрывными. Пространство является пространством скалярных элементов Лагранжа второго порядка. Степени свободы таких элементов связаны со значениями функции в вершинах тетраэдров и в серединах ребер.

Степени свободы векторных элементов имеют смысл для функций и £ Нв(гоЬ\ К), 5 £ (1/2; 1]. В этом случае можно однозначно определить интерполянт яки на каждом тетраэдре К £ Тп. Основной локальный интерполяционный результат приводится в следующей лемме. Лемма 3.3. Пусть и £ Н5(гоС;/С) и 5 > 1/2, имеют место следующие оценки ошибки интерполяции:

||и - ттки\\ок < С/г^||и||лго№ \\rot (и - пки)\\0,к < Ск[{\го1и\5К.

Вводится множество функций, удовлетворяющих дискретному аналогу калибровочного условия

Данное множество является неконформной аппроксимацией множества У0, поскольку <3 \ ^ 0.

В §2 формулируется конечномерная регуляризированная задача. Аналогично непрерывному случаю вводятся операторы £ £(У, и С11 £ как

{В^, Яп) = Ь(У, </„), Уи £ V, Vqгл £

{Снрп, Чн) = с(рл, (?л), Урл, <7Л £ (2л. Определяется билинейная форма а^(-,■)' V х^-» К с параметром /3 > О следующим образом:

1

= а{и,у)+-{Впу,Сн ХВпи) £ V.

Формулируется дискретный аналог задачи (13). Требуется найти такую функцию £ что выполняется равенство:

а$(ил,рл) = - (В^.С^сИУь Г) Уул 6 (17)

Вопросы разрешимости задачи (17) решаются с использованием следующего варианта дискретного разложения Гельмгольца. Лемма 3.4. Для каждой функции V,, £ справедливо разложение:

V* = и^ + Ур„. (18)

Где £ ^о, р^ £ (¡ь- Кроме того, верна оценка:

НЫ^п £ фтог и^Н^ + ||УРй||^п), (19)

где константа С4 не зависит от параметра Л и функции Утверждение о разрешимости дискретной задачи (17) формулируется в следующей теореме.

Теорема 3.1. Решение задачи (17) иЛ £ существует и единственно. Для решения справедлива оценка:

1Ыи,п < С5\\Р\\г

где константа С5 > 0 не зависит от параметра к и Р. В §3 устанавливается оценка сходимости решения дискретной регуляри-зированной задачи (17) к решению обобщенной задачи (13) для случая, когда область П является объединением конечного числа многогранников П^ Считаем, что границы раздела сред дй^ согласованы с разбиением Тк. Теорема 3.2. Пусть м - решение задачи (13) и р - решение задачи (10) обладают дополнительной регулярностью в подобластях Пг: м|Я( £ Н5ЧгоС;/2£), > 1/2, и р|П; £ Я1+5*(П,). <?2 > 1/2. Пусть <ИV £ £2(Ф> а функция ик £ является решением

дискретной регуляризированной задачи (17). В этом случае для ошибки решения верна оценка:

II" - "Jlrot.fi ^ ^НиЦ^вья, +

[ I

где константы С7 > О и С8 > 0 не зависит от параметра к.

В четвертой главе рассматриваются вопросы численного решения дискретной регуляризированной магнитостатической задачи (17). В §1 формулируется система линейных алгебраических уравнений (СЛАУ). Вводятся базисы пространств и (¿ь так, что:

Уц = Брап^о £ = 1,..., ЛГ^} <?л = ярап ¿ = 1 ,...,N<2}.

Формулы 17Л = ^¡^у^!! с//, = Чк<Рк, имеющие место для произвольных функций 6 и € (¿ь, реализуют изоморфизм между пространствами и евклидовыми пространствами вещественных векторов Е^", соответственно. Конечноэлементная обобщенная задача (17) записывается в виде СЛАУ относительно неизвестных коэффициентов разложения и,- как векторов из пространства К"1':

(л+^т)и = Р (20)

Здесь использованы следующие обозначения для векторов и^еЕ"р и матриц А,Т Е X М*":

и = {ио 1 = 1.....щу,

Р = {/¡; к = (Р.Мд - <ВД, Сь-Чпь Р), I = 1.....щ }г,

А = {а,7; аи = а(ЛГ,-, Nj), I,} = 1,..., Ыу }, Т = Ьи = <ВД, С^ВьЦ), I,] = 1.....^ }.

Устанавливается, что Т = ВС~1ВТ, где матрицы В бЕ^х Е^« и С е х определяются как:

В = {Ьу; Ьи = Ь(ЛГ(, <рД I = 1,..., ЛГу, 7 = 1.....л/р },

С = си = с(<р1,<р,), 1,) = 1, ...,N5}.

Свойства матрицы (.А + -^г) формулируются в следующей теореме.

Теорема 4.1. Матрица (л + ^т) симметрична и положительно определена.

Матрицы (а + ^Т^ и (а + ~М^ эквивалентны по спектру

с константами эквивалентности а2 > > 0, не зависящими от параметра Л. Здесь М - матрица масс, соответствующая базису пространства Ук:

М = {ту; ти = ЛГу)о 1,} = 1,..., Щ }.

В §2 предлагается итерационный метод решения задачи (20), скорость сходимости которого не зависит от величины шага сетки. Теорема 4.1. гарантирует, что для решения задачи (20) можно использовать метод сопряженных градиентов с различными вариантами предобуславливания. Строится двухуровневый итерационный метод. На внешних итерациях обращается матрица

^Л + ^ 7^. В качестве предобуславливателя предлагается РеМ = ^Л + ^ М

В свою очередь, для обращения РехС используется метод сопряженных градиентов (внутренние итерации) с ББСЖ предобуславливателем:

(21)

2 - ш \о) ) \о) )

Здесь используется представление (л + ^Д/) = О - £ - и.

Данный метод применялся для решения модельной задачи. Все численные эксперименты проводились на последовательно сгущающихся сетках. Шаги грубой первоначальной сетки сгущались в 2, 3, 4 и в 5 раз. Относительная точность решения внешней векторной задачи составляла гех1 = Ю-10. В то же самое время, точность решения внутренней векторной задачи была существенно снижена и составляла £1п1 = Ю-3. Статистика приводится в следующей таблице.

Таблица 4.

Число рёбер сетки 1364 9672 31364 72880 140660

Время решения, С < 1 4 22 77 198

Число внешних итераций 6 6 6 6 6

Число внутренних итераций 70 200 350 500 750

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

ренних задач В таблице 5 приведены времена решения задач в секундах при различных значениях параметров /? и е1п1 на различных сетках. Таблица 5.

Число рёбер сетки 1364 9672 31364 72880 140660

ß = 10, elnt = Ю-4 < 1 4 25 90 247

ß = l, еш = Ю-3 < 1 4 22 77 198

ß = 0.1, £int = 10"2 1 3 17 55 142

ß = 0.01, е1ш = 10"2 < 1 4 16 48 124

Приведенные в таблице 5 результаты экспериментов показывают, что при уменьшении параметра /? улучшаются спектральные свойства матрицы

(а что в свою очередь приводит к уменьшению количества внутрен-

них итераций. С другой стороны, уменьшение параметра ß приводит к увеличению константы а2 из теоремы 4.1. и как следствие, к увеличению числа внешних итераций. Данные таблицы 5 показывают, что предложенный итерационный метод можно и нужно оптимизировать по параметру ß. Так для параметров ß = 0.01, eint — 10~2 на пятой сетке (число ребер 140660) достигнуто ускорение по сравнению с вариантом из таблицы 4 (/? = 1, eint = Ю-3) более чем в полтора раза. Дальнейшее уменьшение параметра регуляризации приводит к увеличению числа внешних итераций и росту общего времени решения задачи.

В §3 устанавливается возможность использования итерационного метода, предложенного в предыдущем параграфе для решения задач в сложной проводящей среде, содержащей подобласти с большими контрастами по проводимости. В качестве такого примера рассматривается расчетная область, содержащая осадочные породы, зону разлома, соляной диапир и насыщенные нефтью песчаники. Диапир - это конструкция, состоящая из купола и ствола, которая по строению напоминает гриб. Размеры купола составляют 2 ООО х 2 ООО м. в горизонтальных направлениях. Мощность купола (высота) -300 м. Ствол диапира имеет высоту 1200 м. и поперечные размеры 400 х 400 метров. По своим электромагнитным свойствам диапиры характеризуются аномально высокими удельными сопротивлениями. В нашем примере это сопротивление составляло 10 000 Ом м. Расчет такой конструкции интересен с практической точки зрения, поскольку купол диапира образует естественную преграду на пути миграции флюидов (углеводородов) и является

углеводородной ловушкой. Мы учли данную особенность и под куполом разместили насыщенные нефтью песчаники. Суммарная мощность песчаных слоев составила 70 метров, сопротивления слоев варьировалось от 5 до 10 Ом м. Таким образом, контраст по проводимости равен 2 ООО. В нашем примере сетка содержит 197 ООО узлов и 1 150 ООО ребер. Размерность задачи составляет 1 347 ООО скалярных и 2 300 000 векторных неизвестных, относительная точность решения sint — Ю-10. Внешний итерационный процесс сошелся за 19 итераций. Таким образом, можно сделать вывод о том, что предложенный нами итерационный метод позволяет устойчиво получать решение в сложных областях с большими контрастами проводимости.

В заключении приводятся основные результаты, полученные в диссертации.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ

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

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

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

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

ПУБЛИКАЦИИ ПО ТЕМЕ ДИССЕРТАЦИИ

1. Иванов М.И., Катешов В.А., Кремер И.А., Урев М.В. Решение трехмерных стационарных задач импульсной электроразведки //Автометрия. — 2007. - Том. 43. - № 2. - С. 22 - 32.

2. Иванов М.И., Катешов В.А., Кремер И.А., Урев М.В. Решение трехмерных нестационарных задач импульсной электроразведки //Автометрия. -2007. - Том. 43. - № 2. - С. 33 - 44.

3. Иванов М.И., Кремер И.А. Сравнение различных подходов к решению нестационарной системы Максвелла в неоднородной проводящей среде //Тез. Докл. Международная конференция по вычислительной математике 1ССМ - 2009,23 - 25 июня, г. Новосибирск.

4. Иванов М.И., Кремер И.А., Катешов В.А., Эпов М.И. Трехмерное моделирование нестационарных электромагнитных полей от гальванических и индукционных источников: алгоритмы и применение в геоэлектрике //Тез. Докл. 1-ой международной конференции «Актуальные проблемы электромагнитных зондирующих систем» - 2009,27-30 сентября, г. Киев.

5. Иванов М.И., Катешов В.А., Кремер И.А., Эпов М.И. Программное обеспечение модем ЗБ для интерпретации данных нестационарных зондирований с учетом эффектов вызванной поляризации //Записки Горного института - 2009. - Т. 183. - С. 242 - 245.

6. Кремер И.А., Урев М.В. Метод регуляризации стационарной системы Максвелла в неоднородной проводящей среде //Сиб. журн. вычисл. математики. - 2009. - Т. 12,-№2.-С. 161-170.

7. Кремер И.А., Урев М.В. Решение методом конечных элементов регуляри-зированной задачи для стационарного магнитного поля в неоднородной проводящей среде //Сиб. журн. вычисл. математики. - 2010, Т. 13, №1. -С. 33 -49.

Подписано в печать 01.07.2010г. Усл. печ. л. 1,0 Заказ № 125

Формат 60x84 1\16 Тираж 100 экз. •

Отпечатано в типографии ООО « Омега Принт» 630090, г. Новосибирск, пр. Ак.Лаврентьева,6, оф.3-021 тел/факс (383) 335-65-23 email: omegap@yandex.ru

ГЛАВА 1. Введение.

ГЛАВА 2. Обобщенная постановка регуляризированной задачи.

§ 1. Метод регуляризации для абстрактной задачи с седловой точкой.

§ 2. Обобщенная постановка стационарной магнитной задачи.

§ 3. Метод регуляризации для магнитостатической задачи.

ГЛАВА 3. Дискретная регуляризированная магнитостатическая задача

§ 1. Интерполяционные свойства векторных элементов Неделека.

§ 2. Дискретная постановка регуляризированной магнитостатической задачи.

§ 3. Теорема сходимости.

ГЛАВА 4. Численное решение стационарной магнитной задачи.

§ 1. Свойства регуляризированной СЛАУ.

§ 2. Итерационный метод решения регуляризированной стационарной магнитной задачи.

§ 3. Вычисление магнитного поля в проводящей среде с соляным диапиром.

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

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

Рассмотрим систему уравнений Максвелла в квазистационарной постановке, которая часто используется в геофизических приложениях [2]: го£Я=/, сШ?/ = 0, (1.1) го*:Я + |-Д = 0, (ИуВ = 0. (1.2)

Здесь Я и В - вектора напряженности и индукции магнитного поля, Е -напряженность электрического поля, ] - плотность тока. Неизвестные величины связывают уравнения состояния: = аЕ+р, В = 11Н. (1.3)

Здесь - плотность стороннего тока (источник электромагнитного поля), а - удельная электрическая проводимость, [1 — магнитная проницаемость. Предполагаем, что функции а и ц являются скалярными, постоянными в подобластях. На границах раздела подобластей с различными значениями электрической проводимости и магнитной проницаемости выполнены условия сопряжения касательных компонент векторов Я и Е, а также нормальных компонент векторов В и/: тг X Я] = 0, [пх Е] = 0, [В-п] = 0, [/ • п] = 0. (1.4) Здесь п - вектор нормали к границе раздела сред, а квадратные скобки обозначают скачок соответствующей величины при переходе через границу. На различных частях внешней границы, удаленной от источников полей, предполагаем выполнение однородных краевых условий: пхЕ = 0, / • п = 0, п х Я = 0, Вп = 0. (1.5)

Система уравнений (1.1) — (1.3) дополняется начальными данными:

Е\с=0 = Е0, Я|С=0=Я0 или /|г=о=/о> В\^0 = В0. (1.6)

В начальный момент времени поля удовлетворяют следующим стационарным уравнениям:

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

Проведем анализ возможных вариантов решения задачи (1.1) -(1.8). Какие либо варианты симметрий областей не предполагаются, поэтому мы не рассматриваем постановки, которые сводятся к уравнениям меньшей, чем ЗБ, размерности. Известны смешанные постановки задач [3, 4]. в которых в качестве искомых функций выступают одновременно напряженность магнитного поля Н и индукция В. Решение получается в результате минимизации квадратичного функционала основанного на уравнениях состояния (1.3). К недостаткам смешанных постановок можно отнести большое количество неизвестных функций, и как следствие, большие размерности получающихся систем конечномерных уравнений.

Поскольку критической величиной при численном решении трехмерных задач является число неизвестных, одним из возможных подходов к решению системы уравнений (1.1) - (1.8) является исключение одной из функций Е или Я и переход к системе уравнений второго порядка. В случае если исключить напряженность магнитного поля Н, то получим, так называемую, Е - постановку: гоЬ Н0 =/0, rot Е о = О,

Ну /о = О

1.7)

1.8)

17 Вп = 0. о

Такая постановка рассматривалась в работах [5 — 7].

1.9)

При исключении электрического поля Е, ив предположении, что удельная электрическая проводимость а не обращается в ноль, получим Н - постановку [8-13]:

Г 1 дН Js rot — rot Н + — rot~' div 11Н = О, (1Л0)

H\t=0 = Н0.

Уравнение второго порядка получается также, если систему уравнений (1.1) - (1.2) переписать в терминах векторного магнитного потенциала А [14 - 21], [45, 46]:

J.H = rot А.

В этом случае электрическое поле будет выражаться следующим образом: дА

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

Кулона используется в виде div А = 0. В этом случае получается система уравнений, которая связывает неизвестные функции AkU: 1 дА rot — rotA + СТ — + aVU = Js, 11 ot divA = 0, (1.11) дА \ -div a i — + Vt/J + div Js = 0.

Связь А и U в одной системе уравнений вносит дополнительные трудности в решение задачи, которые легко преодолеваются, если калибровочное условие Кулона изменить на следующее условие «Кулоновского типа» [22, 23]: div <тА — 0. (1.12)

Данное изменение калибровочного условия позволяет получить отдельную задачу для скалярного электрического потенциала U. Действительно: дА \ д div стЕ - —div er -—- + VU = — — div а A — div crVU = —div erVU. \ot ) dt

Действуя на первое уравнение системы (1.11) оператором div, выделяем задачу для неизвестной скалярной функции U: div aVU = div Js. (1.13)

Далее считаем функцию U заданной. Получаем следующую задачу для векторного магнитного потенциала с известной функцией U в правой части:

1 дА rot - rot А + er— = -о4U + Is, п л л\ dt div а А = 0.

Выбор одной из систем уравнений (1.9), (1.10) или (1.14) зависит от конкретной постановки задачи. В каждой из систем уравнений имеются свои «минусы», перечислим некоторые из них. Например, в системе (1.9) в правую часть первого уравнения входит производная плотности тока источника по времени djs/dt, а плотность тока часто имеет вид мгновенно включающихся или выключающихся импульсов. В системе (1.10) значение проводимости может быть близко к нулю для плохо проводящих сред, таких как воздух или изоляторы. Прежде чем решать систему уравнений (1.14), необходимо вычислить скалярную функцию U. Тем не менее, все перечисленные системы уравнений имеют одинаковый вид: дХ rot a rot X + ß — = F, (1Л5) div ßX = 0.

Пусть неизвестной функцией будет векторный потенциал А с калибровочным условием (1.12). В пользу данного выбора приведем несколько аргументов. Во-первых, система уравнений (1.14) не содержит множителя 1/<т и слагаемого в правой части. Во-вторых, электрическое поле Е и векторный потенциал А могут иметь один и тот же порядок точности по пространственным переменным, при условии более точного решения задачи относительно неизвестной скалярной функции и. Данное условие легко выполняется, поскольку решение скалярной задачи для и существенно проще решения векторной задачи. В-третьих, нет потери точности и при вычислении измеримых характеристик магнитного поля. Одной из таких характеристик является изменение потока магнитного поля Ф через замкнутый контур Ь. Если поверхность ограниченную контуром Ь обозначить как Б, то формула для вычисления потока выглядит следующим образом: Ф = / ¡лН • п с13. Здесь п - вектор нормали к поверхности 5. Используя формулу Стокса, удается избежать дифференцирования векторного потенциала по пространственным переменным и получить значение потока магнитного поля без потери точности. Действительно:

Таким образом, дальнейшие исследования проводятся для векторного потенциала А, но все полученные результаты без ограничения переносятся и на случаи Е и Н постановок (1.9), (1.10).

Определим объект наших исследований. В начальный момент времени векторный потенциал удовлетворяет следующей системе уравнений: 1 rot — rotA fi

UL Лд — —U V div аА0 = 0,

-oW + Js,

1.16)

Решение магнитостатической задачи осложняется тем, что оператор rot обладает нетривиальным ядром, которое содержит градиенты скалярных функций. Единственность решения обеспечивает калибровочное условие. Кроме этого, необходимо выполнение условия div{—aVU + Js) = О, являющегося необходимым для разрешимости задачи (1.16). Данное условие также называется условием совместности правой части. Заметим, что в нашем случае выполнение условия совместности гарантировано в силу определения скалярного потенциала U как решения задачи (1.13).

Сравним системы уравнений (1.14)и(1.16)с точки зрения методов их решения. Введем временную сетку и перейдем от задачи (1.14) тем или иным способом к дискретной по времени постановке [24 - 26]. На одном временном шаге общий вид задачи будет следующим: 1 rot — rot А + ааА = F, ' п л

1 К1-1-') div о А = 0.

Число а > 0 зависит от величины шага по времени и от выбранной схемы интегрирования (для задач в частотной области зависимость от гармонической частоты), а правая часть системы F дополнительно будет зависеть и от значений векторного потенциала А, вычисленных на соседних временных шагах. При переходе к обобщенным постановкам в пространстве H(rot; Д) билинейная форма, соответствующая левой части первого уравнения системы (1.17) эллиптична, слагаемое ааА обладает регуляризирующим свойством. Данное обстоятельство позволяет вместо системы (1.17) решать только первое уравнение. Выполнение калибровочного условия для векторного потенциала А является следствием соленоидальности правой части div F = 0. Для стационарной задачи (1.16) аналогичные рассуждения в пространстве H(rot; П) не работают. При переходе к обобщенной постановке билинейная форма, порожденная левой частью первого уравнения системы (1.16), не является положительно определенной, а только неотрицательной. В этом случае стационарную . задачу можно рассматривать только как систему уравнений (нельзя игнорировать калибровочное условие), и в этом- смысле стационарная задача (1.16) является более сложной, чем нестационарная задача (1.17). Косвенно данную мысль подтверждает и сравнение количества публикаций по решению стационарных и нестационарных (во временной и частотной области) задач. Разница на порядки в пользу нестационарных задач.

С другой стороны, отсутствие соленоидальности правой части F первого уравнения системы (1.17) с необходимостью ведет к нарушению калибровочного условия для векторного потенциала А. Поскольку F определяется численными методами, то (Ну ¥ Ф 0 получается в результате накопления погрешностей вычислений. Особенно сильно данные эффекты проявляются при решении задач на установление полей (методы переходных процессов, электромагнитные зондирования становлением полей) на поздних временах и разрушают численное решение [27, 28]. Если задача содержит подобласти с большими контрастами по проводимости, то неустойчивости решений могут проявляться и на ранних временах. Отметим также, что аналогичная ситуация возникает и при решении прямых задач в частотной области. При малых значениях гармонической частоты ухудшаются регуляризирующие свойства слагаемого ааА, а соответственно растет число обусловленности системы линейных алгебраических уравнений [29].

В предельном варианте а = 0 все перечисленные проблемы возникают при решении стационарной задачи (1.16). Очевидно, что если будет решена стационарная задача, то будет понятно каким образом преодолевать проблемы и в нестационарном случае. Актуальность представленной работы обоснована.

Объектом наших исследований будет стационарная система уравнений (1.16). Нам потребуются постановки, допускающие несоленоидальные правые части div F Ф 0, и гарантирующие выполнение калибровочного условия (1.12). Предметом исследований данной работы будет решение таких систем методом векторных конечных элементов.

Анализ источников и степень разработанности проблемы. Из литературы известно несколько способов удовлетворения калибровочному условию. Наиболее простой из них - ничего специального не предпринимать. Так в работах [47 - 49] отмечено, что можно рассматривать решение вырожденного rot — rot уравнения методом сопряженных градиентов. Замечено, что одна итерация данного метода не портит калибровочного условия. Иными словами, если начальное приближение удовлетворяет, то и все последующие приближения будут удовлетворять ограничениям на дивергенцию от решения. Существенным недостатком такого подхода является то, что вычисления должны быть абсолютно точными. И вообще говоря, непонятно, почему должен сходиться метод сопряженных градиентов.

Калибровочное условие можно рассматривать как ограничение на решение. В этом случае данное ограничение может быть учтено введением в постановку задачи дополнительной скалярной неизвестной функции р - множителя Лагранжа. Одной из первых работ на эту тему является теоретическая работа [34], в которой рассматривается полная система уравнений Максвелла в частотной области. Постановки задач с множителем Лагранжа допускают использование несовместных правых частей. В работе [24] метод множителей Лагранжа применялся для решения стационарной и нестационарной систем уравнений Максвелла во временной области* с разрывными коэффициентами. В работе [39] рассмотрены постановки для нестационарной системы уравнений. В публикации [53] постановка с множителем Лагранжа приводится для задачи на собственные значения. Такие постановки приводят к задаче с седловой точкой. В операторном виде задача записывается в виде системы уравнений для неизвестных функций (Ар):

1,1) (сЛ), diva - блок (2,1) (В), crV - блок (1,2) (Ът). Общая теория решения задач с седловой точкой приводится, например, в работах [29, 30]. Для численного решения систем линейных алгебраических уравнений, соответствующих задачам с седловой точкой, зачастую применяются методы, позволяющие выделить отдельные задачи для неизвестных функций и в которых требуется вычислять дополнение Шура [54]. В случае стационарной системы уравнений Максвелла применение таких методов осложнено тем, что диагональные блоки (1,1) и (2,2) вырождены. В качестве альтернативного подхода можно решать полную систему уравнений, не производя разделение неизвестных функций. Так же можно применить прием регуляризации [29], суть которого заключается в том, что блок (2,2) меняется на невырожденный —£ V, но содержащий малый параметр регуляризации - £. После таких изменений исходной матрицы удается разделить переменные и заменить одну большую задачу на две меньшей размерности. Минус подхода заключается в том, что параметр регуляризации должен быть мал, что ведет к неустойчивости вычислений, поскольку блок (2,2) становится «почти вырожденным». Общим плюсом постановки с множителем Лагранжа является то, что калибровочное условие в явном виде

1.18) 1

Здесь оператор rot- rot в исходной задаче (1.16) соответствует блоку включается в систему уравнений и полученное решение, естественно, удовлетворяет ему вне зависимости от соленоидальности правой части.

Следующий способ получения невырожденной системы уравнений основан на известной формуле векторного анализа и называется "grad — div" регуляризацией. Так в однородном случае {¡х = а = 1) имеем rot rot и — V(div и) = —V2u. (1-19)

В правой части данного тождества находится положительноопределенный оператор (при подходящей области определения оператора и соответствующих краевых условиях). В наиболее простом [57] варианте для задачи о векторном потенциале и, с учетом того, что div и = 0, оператор rot rot заменяется оператором —V2. Далее без проблем решается эллиптическая задача. Существенным ограничением такого подхода является то, что он может применяться только для однородных электромагнитных сред. Но самое главное заключается в том, что после такой замены операторов решение не обязано удовлетворять калибровочному условию. Магнитное поле при этом восстанавливается через векторный потенциал однозначно. В другом варианте [58] к оператору rot rot добавляется регуляризирующее слагаемое —V(div и). Таким образом, снова получается задача с положительноопределенным оператором. Минусы такого подхода выявляются при более строгом анализе регуляризированной задачи. Фактически, в этом случае возникают дополнительные требования на решение: требуется непрерывность не только касательных компонент, но также и нормальных. Эти требования позволяют определять только непрерывные решения.

Для того чтобы преодолеть ограничения однородности, следует расширить понятие дивергенции от решения, полагая, что div и Е Y(£i) i i для некоторого подходящего Гильбертова пространства, К(П): Такие расширения приводят к обобщенным постановкам с билинейной формой rot и, rot v)x + ос (div и, div v)y. (1.20)

Здесь скобки означают скалярное произведение в Гильбертовых пространствах X(Q) (как правило, X = L2 (&■)) и Y(й) соответственно, а > 0 некоторый параметр регуляризации. Примеры таких расширений приведены в работе [59] и в серии работ [55, 56]. Следует помнить, что расширение понятия дивергенции автоматически влечет необходимость адекватной аппроксимации пространства Y и формы (div и, div v)Y • Более полную библиографию по вопросам "grad — div" регуляризаций можно найти в обзорной работе [50].

Еще один вариант учета — это разрывный метод Галеркина (DG) или метод внутренних штрафов (interior penalty). Дело в том, что после дискретизации области на элементы, на гранях этих элементов можно определить поверхностную дивергенцию через скачки нормальных компонент функций. Далее в билинейной форме (1.20) вместо скалярного произведения (div и, div v)Y> определенного в объеме, можно использовать его поверхностный аналог. Теория DG для решения системы уравнений Максвелла приводится в работах [51, 52]. Получается, что на каждой внутренней грани разбиения области необходимо иметь два значения функции: «справа» и «слева» от грани. Данное обстоятельство автоматически приводит к увеличению размерности системы линейных уравнений, со всеми вытекающими последствиями. Дополнительно следует исследовать вопрос о величине параметра регуляризации а в формуле (1.20).

Следующая серия работ посвящена теме дискретной регуляризации. Суть метода дискретной регуляризации заключается в способах численного решения конечномерной задачи, соответствующей системе уравнений (1.18). Можно выделить два способа такой регуляризации. Введем обозначения для конечномерных операторов и. неизвестных функций <Ah, 3h, Dh, Ah, ph. Первый способ (weighted augmentation) основывается на том, что ЪкАк = 0. В этом случае для некоторой весовой матрицы Wh можно перейти к системе уравнений

Если при этом оказывается, что матрица (<Ah + положительно определена, то к решению данной системы уравнений можно применять различные способы, основанные на невырожденности блока (1,1). В работах [63, 64] предлагаются методы решения (1.21), основанные на изучении алгебраических свойств дискретных линейных операторов <Ah, Ън, Wh . Отметим, что система (1.21) по своей сути остается задачей с седловой точкой, решение которой требует многократное обращение блока (1,1). Второй способ решения дискретного аналога задачи (1.18) основан на регуляризации блока (2,2) и решении приближенной задачи [50]

Ah Ъ/ \ (AhE\ = (Fh\ -EVhJ\pheJ У or

В этом случае имеется возможность произвести исключение неизвестной функции phe и перейти к задаче меньшей размерности для функции Ah£

В работе [63] такой способ был реализован для s = 1 и T>h -конечномерного аналога оператора div(V). Выбранные значения параметров s и T)h имеют смысл, если ph = 0 (что верно в случае div Fh = 0). Несомненное достоинство такого подхода заключается в том, что матрица обращается только один раз.

1 —1

Заметим, в случае

Wh = -Vh~x блок (1.1) задачи (1.21) совпадает с матрицей из второго способа дискретной регуляризации. Данный способ допускает обобщение на случай div F Ф 0, а также может применяться для неоднородных по электромагнитным свойствам сред.

Рассмотрим теоретические основы используемых далее подходов. При формулировке обобщенных постановок, наиболее естественным представляется сосредоточиться на использовании векторного функционального пространства H(rot; П), в котором основным свойством гладкости элементов является непрерывность их касательных компонент, а данное свойство является фундаментальным свойством искомого решения задачи (1.16). Отметим несколько работ, освещающих различные аспекты такого использования пространства Я (rot; £>) и полезных для наших дальнейших исследований. В работе [60] приводится точная характеризация пространства следов HQrot; £2) и предлагается конструктивный метод построения непрерывного оператора продолжения. В работах [35, 65] при дополнительных ограничениях на дивергенцию устанавливаются критерии компактности ограниченных множеств функций в пространстве L2(fl). При исследовании вопросов разрешимости обобщенных задач важным инструментом являются различные варианты разложения Гельмгольца и соответствующие оценки для таких разложений. Отметим работы [29, 61, 65], в которых выполнены разложения и получены оценки для однородных случаев.

Ключевыми вопросами численного моделирования являются вопросы аппроксимации пространств и теория интерполяции. Отметим работы Неделека [36, 62], в которых описываются векторные конечные элементы, порождающие конформные аппроксимации пространства HQrot; Q) и приводятся основные факты теории интерполяции для функционального пространства H(rot; Q).

Цели данной работы следующие.

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

• Методом векторных конечных элементов Неделека получить конечномерную постановку регуляризированной задачи и исследовать ее.

• Установить сходимость конечноэлементного решения к точному решению обобщенной регуляризированной задачи.

• Получить систему линейных алгебраических уравнений (СЛАУ), установить ее свойства. Предложить и обосновать численный метод решения СЛАУ.

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

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

Теоретическая значимость представленной работы заключается в постановке эквивалентной регуляризированной задачи для стационарного магнитного поля в неоднородной проводящей среде. Постановка допускает несовместные правые части. Калибровочное условие является свойством решения регуляризированной задачи. Для решения регуляризированной задачи применяется метод векторных конечных элементов. Доказана теорема сходимости конечноэлементного решения к точному. Предложенная регуляризированная постановка задачи позволяет получить СЛАУ с симметричной и положительно определенной матрицей.

Прикладная ценность работы заключается в том, что регуляризированные постановки задач технологично решаются векторным методом конечных элементов. Полученные результаты распространяются на нестационарные постановки (как в частотной области, так и во временной). Результаты исследований реализованы в виде численных модулей и включены в программу «MODEM 3D» [41], предназначенную для интерпретации данных 3D нестационарных зондирований с учетом эффектов вызванной поляризации.

Личный вклад автора. Основные результаты работы получены автором. Достаточные условия теоремы 2.3. установлены совместно с М.В. Уревым. Постановку практической задачи, приведенной в параграфе 3 четвертой главы, сделал М.И. Эпов. Программная реализация «MODEM 3D» выполнена совместно с М.И. Ивановым. Представленные и выносимые на защиту результаты, полученные в совместных исследованиях, согласованы с соавторами.

Апробация работы. Результаты работы докладывались на семинарах ИВМ и МГ СО РАН (2008, 2009, 2010, г. Новосибирск), Baker Atlas (2008, г. Новосибирск). Также на Российских и международных семинарах и конференциях.

6-й Международный геофизический научно-практический семинар Применение современных электроразведочных технологий при поисках месторождений полезных ископаемых, 25-27 марта 2008 г., г. Санкт - Петербург.

5-ая Международная конференция и выставка «Недра - 2008», 1-4 апреля 2008 г., г. Москва.

Научно - практическая конференция Комплексирование геолого - геофизических методов при обосновании нефтегазопоисковых объектов на Сибирской платформе (в Восточной

Сибири и Республике Саха (Якутия)), 21-23 апреля 2008 г., г. Новосибирск.

Международная конференция «Математические Методы в Геофизике - 2008», 13-15 октября 2008 г., г. Новосибирск.

2-ая Всероссийская конференция «Актуальные проблемы строительной отрасли», 14-15 апреля 2009 г., г. Новосибирск.

Международная конференция по вычислительной математике 1ССМ -2009, 23 - 25 июня, г. Новосибирск.

1-ая международная конференция «Актуальные проблемы электромагнитных зондирующих систем», 27 - 30 сентября 2009 г., г. Киев.

Публикации. Основные результаты диссертации опубликованы в 7 работах, в том числе 5 работ в журналах из перечня ВАК [22, 23, 41, 42,

43] и две работы [27, 28] - тезисы международных конференций.

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

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

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

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

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

Структура и логика дальнейшего изложения следующая. Работа, помимо введения, содержит еще три главы, заключение и список литературы. Во второй главе в первом параграфе приводится схема решения абстрактной задачи с седловой точкой. Формулируются достаточные условия ее разрешимости. Рассматривается прием регуляризации, который позволяет разделить неизвестные функции и получить отдельные приближенные формулировки задач для каждой из искомых функций. Данные результаты имеют справочный характер и позволяют точно сформулировать проблемы, возникающие при решении задач с седловой точкой. Во втором параграфе стационарная система уравнений для векторного потенциала (1.16) формулируется как обобщенная задача с седловой точкой. В качестве неизвестных функций выступают векторный потенциал и множитель Лагранжа. Производится проверка выполнения условий разрешимости обобщенной задачи. Для этих целей исследуется структура подпространства функций, удовлетворяющих калибровочному условию Кулоновского типа (1.12).

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

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

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

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

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

ГЛАВА 2. Обобщенная постановка регуляризированной задачи

Выводы по главе 4. Сформулирована система линейных алгебраических уравнений, соответствующая дискретной регуляризированной задаче. Исследованы свойства матрицы СЛАУ: установлена ее симметричность, положительная определенность и внутренних итераций необходимых для обращения матрицы спектральная эквивалентность более простой матрице (а + ^М^. Для метода сопряженных градиентов предложен оптимальный предобуславливатель. Оптимальность заключается в том, что число итераций метода не зависит от шага сетки. На примере численного решения модельной задачи установлено, что параметр регуляризации /? может быть использован для оптимизации предложенного метода. Показано, что скорость решения регуляризированной стационарной задачи сравнима со скоростью расчета одного временного шага нестационарной задачи без регуляризации. Данный факт и теорема 4.1. являются основными результатами, полученными в четвертой главе. С помощью предложенного итерационного метода решена практическая задача в трехмерной области со сложной геометрией и большими контрастами по проводимости.

ЗАКЛЮЧЕНИЕ

Сформулируем основные результаты, полученные в диссертации.

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

2. Предложен способ регуляризации трехмерной стационарной магнитной задачи в неоднородной проводящей среде. Сформулирована и исследована обобщенная регуляризированная задача для векторного магнитного потенциала.

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

4. Обобщенная регуляризированная задача сформулирована и исследована в конечномерных пространствах Неделека на тетраэдральных конечноэлементных сетках.

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

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

7. Предложен оптимальный итерационный метод решения СЛАУ, скорость сходимости которого не зависит от величины шага сетки.

Эффективность метода установлена на примере численного решения модельной задачи.

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

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

Автор выражает глубокую признательность своему научному руководителю М.В. Уреву за внимание и помощь в работе. Автор также благодарен Е.Ю. Антонову, М.И. Иванову, В.П. Ильину, В.А. Катешову, B.C. Могилатову, М.И. Эпову за полезные обсуждения и плодотворное сотрудничество.

1. Могилатов B.C., Балашов Б.П. Зондирование вертикальными токами (ЗВТ) //Изв. РАН Сер.: Физика Земли, №6 (1994), С.73.

2. Тамм И.Е. Основы теории электричества, М.: Наука, 1989, 504 с.

3. Bossavit A. Computational Electromagnetism: Variational Formulations, Complementarity, Edge Elements, Academic Press, San Diego, 1998.

4. Albanese R., Rubinacci G. Analysis of three-dimensional electromagnetic fields using edge elements //J. Comput. Phys. 1993. — Vol. 108.-P. 236-245.

5. Alonso A. and Valli A. A domain decomposition approach for heterogeneous time-harmonic Maxwell equations //Comput. Methods Appl. Mech. Engrg. 1997. - Vol. 143. - P. 97-112.

6. Bermudez A., Rodriguez R., and Salgado P. Numerical analysis of an electric formulation of the eddy current problem //Numer.Math. 2005. Vol. 102. P. 181-201.

7. Golias N.A., Anntonopoulus C.S., Tsiboukis T.D., and Kriezis E.E. 3D eddy current computation with edge elements in terms of the electric intensity //COMPE. 1998. - Vol. 17. - P. 667-673.

8. Alonso A., and Valli A. An optimal domain decomposition preconditioner for low-frequency time-harmonic Maxwell equations //Math. Сотр. 1999. - Vol. 68. - P. 607-631.

9. Alonso A., Hiptmair R., and Valli A. Mixed finite element approximation of eddy current problems// IMA J. Numer. Anal. 2004. -Vol. 24-P. 255-271.

10. Bermudez A., Rodriguez R., and Salgado P. A finite element method with Lagrange multipliers for low-frequency harmonic Maxwell equations //SIAM J. Numer. Anal. 2002. - Vol. 40. - P. 1823-1849.

11. Bermudez A., Rodriguez R., and' Salgado P. Numerical solution of eddy current problems in bounded domains using realistic boundary conditions //Comput. Methods Appl. Mech. Engrg. 2005. - Vol. 194. -P. 411-426.

12. Bermudez A., Rodriguez R., and Salgado P. Numerical treatment of realistic boundary conditions for the eddy currents problem in an electrode via Lagrange multipliers //Math. Comp. 2005. - Vol. 74. - P. 123-151.

13. Tu H.T., Shao K.R., and Zhou K.D. H method for solving 3D eddy current problems //IEEE Trans. Magn. 1995. - Vol. 31. - P. 35183520.

14. Albertz D., and Henneberger G. Calculation of 3D eddy current fields using both electric and magnetic vector potential in conducting regions //IEEE Trans. Magn. 1998. - Vol. 34. - P. 2644-2647.

15. Albertz D., and Henneberger G. On the use of the new based A-A,T formulation for the calculation of time-harmonic stationary and transient eddy current field problems //IEEE Trans. Magn. 2000. Vol. 36. - P. 818-822.

16. Biro O. Edge element formulations of eddy current problems //Comput. Meth. Appl. Mech.Eng. 1999. - Vol. 169. - P. 391-405.

17. Biro O., and Preis K. On the use of the magnetic vector potential in the finite element analysis of the three-dimensional eddy currents //IEEE Trans. Magn. 1989. - Vol. 25. - P. 3145-3159.

18. Biro O., and Valli A. The Coulomb gauged vector potential formulation for the eddy-current problem in general geometry: well-posedness and numerical approximation //University di Trento. 2005. - UTM 682 preprint.

19. Leonard P.J., and Rodger D. Finite element scheme for transient 3D eddy currents //IEEE Trans. Magn. 1988. - Vol. 24. - P. 90-93.

20. Morisue T. Magnetic vector potential and electric scalar potential in three-dimensional eddy current problem //IEEE Trans. Magn. 1982. -Vol. 18.-P. 531-535.

21. Morisue T. 3D-eddy current calculation using the magnetic vector potential //IEEE Trans. Magn. 1988. - Vol. 24. - P. 106-109.

22. Иванов М.И., Катешов B.A., Кремер И.А., Урев М.В. Решение трехмерных стационарных задач импульсной электроразведки //Автометрия. 2007. - Том. 43. - № 2. - С. 22 - 32.

23. Иванов М.И., Катешов В.А., Кремер И.А., Урев М.В. Решение трехмерных нестационарных задач импульсной электроразведки //Автометрия. 2007. - Том. 43. - № 2. - С. 33 - 44.

24. Chen Z., Du Q., and Zou J. Finite element methods with matching and nonmatching meshes for Maxwell equations with discontinuous coefficients //SIAM J. Numer. Anal. 2000. - Vol. 37. - P. 1542-1570.

25. Ciarlet P. Jr. and Zou J. Fully discrete finite element approaches for time-dependent Maxwell's equations //Numer. Math. 1999. - Vol. 82. -P. 193-219.

26. Makridakis Ch. G., and Monk P. Time-discrete finite element schemes for Maxwell's equations //M2AN Math. Modelling and Numer. Anal. -1995.-Vol. 29.-P. 171-197.

27. Иванов М.И., Кремер И.А. Сравнение различных подходов к решению нестационарной системы Максвелла в неоднородной проводящей среде //Тез. Докл. Международная конференция по вычислительной математике ICCM — 2009, 23 25 июня, г. Новосибирск.

28. Girault V., Raviart Р.-А. Finite element methods for Navier Stokes Equations. Springer - Verlag, Berlin, 1986.

29. Brezzi F., Fortin M. Mixed and Hybrid Finite Element Methods. Springer Verlag, New-York, 1991.

30. Лионе Ж.-Л., Мадженес Э. Неоднородные граничные задачи и их приложения. М.: Мир, 1971.

31. Дюво Г., Лионе Ж.-Л. Неравенства в механике и физике. — М.: Наука, 1980.

32. Kikuchi F. On a discrete compactness property for the Nedelec finite elements //J. Fac. Sci. Univ. Tokyo Sect. LA Math. 1989. - Vol. 36. -№ 3 - P. 479-490.

33. Гудович И. С., Крейн С. Г., Куликов И. М. Краевые задачи для уравнений Максвелла //Доклады АН СССР. Т. 207. - №2. - С. 321 -324.

34. Weber С. A local compactness theorem for Maxwell's equations //Math. Meth. In the Appl. Sci. 1980. - Vol. 2. - P. 12-25.

35. Nedelec J.C. A new family of mixed finite elements in R3 //Numer. Math. 1986. - Vol. 50. - P. 57-81.

36. Сьярле Ф. Метод конечных элементов для эллиптических задач. М.: Мир, 1982.

37. Amrouche С., Bernard! С., Dauge М. and Girault V. Vector Potentials in Three-dimensional Non-smooth Domains //Math. Meth. Appl. Sci. -1998.-Vol. 21.-P. 823-864.

38. Assous F., Degond P., Heintze E., Raviart P., Segre J. On finite-element method for solving the three-dimensional Maxwell equations //J. Comput. Phys. 1993. - Vol. 109. - P. 222 - 237.

39. Самарский А.А., Лазаров Р.Д., Макаров В.Л. Разностные схемы для дифференциальных уравнений с обобщенными решениями, М.: Высш. Шк., 1987 296 с.

40. Иванов М.И., Катешов В.А., Кремер И.А., Эпов М.И. Программное обеспечение модем 3D для интерпретации данных нестационарных зондирований с учетом эффектов вызванной поляризации //Записки Горного института 2009. - Т. 183. - С. 242 - 245.

41. Кремер И.А., Урев М.В. Метод регуляризации стационарной системы Максвелла в неоднородной проводящей среде //Сиб. журн. вычисл. математики. 2009. - Т. 12. - №2. - С. 161 - 170.

42. Кремер И.А., Урев М.В. Решение методом конечных элементов регуляризированной задачи для стационарного магнитного поля в неоднородной проводящей среде //Сиб. журн. вычисл. математики. -2010, Т. 13, №1.-С. 33-49.

43. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. М.: Мир, - 1991. - 367 с.

44. Albanese R., Rubinacci G. Magnetostatic field computations in term of two-component vector potentials //Int. J. Num. Meth. In Eng. 1990. -vol. 29. - P. 515-532.

45. Biro O., Preis K., Richter K.R. On the use of the magnetic vector potential in the nodal and edge finite element analysis of 3D magnetostatic problems //IEEE. Trans. Magn. 1996. - vol. 32. - P. 651 -654.

46. Ren Z. Influence of the R.H.S. on the convergence behaviour of the curl — curl equation //IEEE. Trans. Magn. 1996. - Vol. 32. - №3. --P. 655-658.

47. Ren Z., Ida N. Derivation of various dual formulations in magnetostatics via error based energy approach //IEEE. Trans. Magn. 1999. - Vol. 35. — №3. - P. 1167-1170.

48. Bouillault F., Buffa A., Maday Y., Rapetti F. The mortar edge element method in three dimensions: application to magnetostatics //SIAM J. Sci. Comput. 2003. - Vol. 24. - №4. - P. 1303 - 1327.

49. Hiptmair R. Finite elements in computational electromagnetism //Acta Numerica. 2002. - Vol. 11. - P. 237 - 339.

50. Perugia I., Schotzau D., Monk P. Stabilized interior penalty methods for the time-harmonic Maxwell equations //Comp. Meth. Appl. Mech. Engrg. Vol. 191. -P. 4675 - 4697.

51. Kikuchi F. Mixed and penalty formulations for finite element analysis of an eigenvalue problem in electromagnetism //Computer Meth. Appl. Mech. Eng. 1987. - Vol. 64. - P. 509 - 521.

52. Kikuchi F., Hara M., Wada T. A mixed finite element method for 3-D analysis of cavity resonators //Sci. Papers. College. Arts. Sci. Univ. Tokyo. 1992. - Vol. 42. - P. 111 - 129.

53. Benzi M., Golub G., Liesen J. Numerical solution of saddle point problems //Acta Numerica. 2005. - Vol. 14. - P. 1 - 137.

54. Bramble J., Pasciak J. A new approximation technique for div-curl systems //Math. Comp. 2004. - Vol. 73. - №248. - P. 1739-1762.

55. Bramble J., Kolev T., Pasciak J. The approximation of the Maxwell eigenvalue problem using a least-squares method //Math. Comp. 2004. -Vol. 74.-№252.-P. 1575-1598

56. Haber E., Ascher U., Aruliah D., Oldenburg D. Fast simulation of 3D electromagnetic using potentials //J. Comput. Phys. 2000. — Vol. 163. -P. 150-171.

57. Haber E., Ascher U. Fast finite volume simulation of 3D electromagnetic problems with highly discontinuous coefficients //SIAM J. Sci. Comput. 2001. - Vol. 22. - №6. - P. 1943-1961.

58. Costabel M., Dauge M. Weighted regularization of Maxwell equations in polyhedral domains //Numer. Math. 2002. - Vol. 93. - P. 239 - 277.61.62. 63. 63.