Mоделирование электромагнитных полей в среде с анизотропной электропроводностью Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 11, № 3, 2006

ЫОДЕЛИРОВАНИЕ ЭЛЕКТРОМАГНИТНЫХ ПОЛЕЙ В СРЕДЕ С АНИЗОТРОПНОЙ ЭЛЕКТРОПРОВОДНОСТЬЮ*

Н.В. Орловская, Э.П. Шурина Новосибирский государственный технический университет, Россия e-mail: [email protected], [email protected]

М. И. Эпов

Институт геофизики СО РАН, Новосибирск, Россия e-mail: [email protected]

Using the method of differential forms we have shown that for modeling of electromagnetic fields in the anisotropic medium the system of the first order Maxwell's equations is preferable compared to the second order equations. Results of the numerical tests are presented.

1. Уравнения Максвелла

Рассмотрим систему уравнений Максвелла в однородной области:

' dtB = -V х Е, dtD = V х Н - аЕ - J, div(D) = р, div(B) = 0.

Здесь е, а — скалярные величины. Уравнения второго порядка для электрического и магнитного полей имеют следующий вид:

, 1 d2D д Е dJ

д 2еН д

rot (rot Н) + -—— + а— (juH) - rot(J) = 0. dt2 dt

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

H(rot;Q) = |v G [L2(tt)]3 : rotv G [L2(^)]3}

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-05-64528), совместного международного проекта NWO и Российского фонда фундаментальных исследований (грант № 047.016.003).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

При аппроксимации системы уравнений Максвелла как системы первого порядка вместе с пространством H(rot, П) вводят пространство H(div, П):

H(div; П) = jv е [Ь2(П)]3 : divv е ¿2(П)} .

Функции из пространства H(div; П) аппроксимируют векторные величины потока, такие как D, B, функции из пространства H(rot; П) — векторные поля E, H. Для пространств H1 (П), H(rot, П), H(div, П) выполняются отношения включения, называемые комплексом де Рама [1]:

H 1(П) —^ H(rot; П) H(div; П) Ь2(П).

Введем обозначения для дискретных подпространств пространств H 1(П), H(го^П), H(div; П):

— пространство непрерывных кусочно-полиномиальных функций

Б(П) = {ф е H 1(П)};

— пространство Неделека роторно-конфорных функций [2, 3]

ТО(П) = {v е H(rot; П)};

— пространство Равьяр — Тома дивергентно-конформных функций [4]

RT(П) = {u е H(div; П)};

— пространство полиномиальных функций

Q^) = {w е Ь2(П)}.

Для дискретных пространств S(П), ЖД(П), RT(П), Q(П) комплекс де Рама принимает вид

S(П) —^ ТО(П) RT(П) Q(П).

Параметры е, ^, о не всегда являются скалярными величинами. В общем случае это положительно определенные тензоры второго ранга. Исследование уравнений Максвелла с тензорными коэффициентами удобнее всего проводить, используя аппарат дифференциальных форм [1]. Дифференциальные формы известны давно и первоначально были предложены в качестве удобного аппарата в физике и дифференциальной геометрии [5]. Для дифференциальных форм в R3 введем следующие обозначения:

0-форма f,

1-форма f1dx + f2dy + f3dz,

2-форма f1dy Л dz + f2dz Л dx + f3 dx Л dy,

3-форма fdx Л dy Л dz.

Дифференциальные формы обладают следующими свойствами [1, 5-7].

Таблица 1. Соотношения между функциональными пространствами и дифференциальными формами

Функциональное Дифференциаль- Внешний Конечно-

пространство ная форма дифференциал элементные

пространства

H1 (Q) 0-форма grad S (Q)

H(rot; Q) 1-форма rot ND(Q)

H(div; Q) 2-форма div RT (Q)

L2(Q) 3-форма 0 Q(Q)

Для 1-формы и выполняется

J ивГ = J (¡, т)вГ,

г г

где Г — гладкая кривая с касательным вектором т. Для 2-формы и имеем

У твБ = J (¡,п)в8,

в 5

где 8 — ориентированная площадь с внешней нормалью п. Определим оператор внешнего дифференциала в [8]

i= 1

dfi(x) j л ^ dfp(x)

- dxi Л uXi-, Л ... Л dxip + ... + / - axi л ^х^-, dxi 1 p ^ dxi

dxj, Л dx^, Л ... Л dxjp,

i= 1

где P — порядок дифференциальной формы.

Обозначим * оператор Hodge [9], действующий следующим образом:

*(dxi1 Л ... Л dxip) = d.xip+1 Л ... Л dxin.

Для внешнего дифференциала выполняются тождества

u = d<ф = grad ф, w = du = rot u, v = dw = div w,

где ф — 0-форма, u — 1-форма, w — 2-форма, v — 3-форма. Для оператора Hodge * выполняются следующие тождества:

*ф = v, *u = w, *w = u, *w = ф,

n

где ф — 0-форма, u — 1-форма, w — 2-форма, v — 3-форма. Оператор Hodge обеспечивает выполнение материальных соотношений, где переменные E, H являются 1-формами, а B, D — 2-формами.

Соотношения между функциональными пространствами и дифференциальными формами представлены в табл. 1.

2. Уравнение второго порядка для электрического поля

В данной работе рассмотрены среды с анизотропными проводящими свойствами, в которых параметр a — тензор ранга 2. Входящие в систему уравнений Максвелла переменные являются следующими дифференциальными формами: E, H — дифференциальные формы первого порядка (1-формы); B, D, J — дифференциальные формы второго порядка (2-формы).

Закон Ампера в дифференциальных формах:

f = -- (1)

Закон Фарадея в дифференциальных формах:

dD

д- = dH - a*E - J. (2)

Материальные соотношения принимают вид

D = e*E; (3)

B = ; (4)

J = a*E. (5)

Будем полагать, что — диагональные тензоры, а a — тензор ранга 2. Умножим уравнение (1) на

дВ

— = -^-1*dE. (6)

Подставив (4) в (6) и выполнив преобразования в соответствии со свойствами оператора Hodge, получим

dH

^ = -^-1*dE. (7)

Поскольку (7) содержит 1-формы, взятие внешнего дифференциала от этого уравнения равносильно взятию операции rot. Вычислим внешний дифференциал от (7):

^ = -d(/'-1*dE).

(8)

Подставив (2) в (8) и преобразовав, получим

д2 D д д J

d("-1*dE) + Ж + 8t (a*E ) = - Ж' (9)

д

Преобразуем — (a*E):

д д

— (a*E) = — [(ajEx + a^Ey + a^Ez)dy Л dz + (a^E + a^y + a^)dz Л dx+

/ д д д \ + (a?Ex + a23Ey + a33Ez)dx Л dy] = i — №*) + ^(a^Ey) + ^(a^Ez ) j dy Л dz+

д д д д д + ( д^ (a2Ex) + ^ (a2Ey) + - (a2Ez)) dz Л dx + (- (a2Ex) + ^ (a2Ey) +

+ Ji(a|EzЛ dx Л dy = дО"*E + a*(10)

Тогда уравнение (9) принимает вид

WdE) + W + Tt *E + а*-dt = - Tt. (11)

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

[n х E]ria = 0; (12)

[n ■ аЕ]Г12 = 0. (13)

Использование векторных функций из пространства H(rot; П) обеспечивает выполнение условия (12). Условие (13) следует учитывать в вариационной постановке в виде дополнительного дивергентного ограничения.

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

3. Уравнение второго порядка для магнитного поля

Для перехода к уравнению второго порядка для магнитного поля в анизотропной среде умножим уравнение (2) на е-1*:

dD

£-1т = £-1*dH - е-1*(а*Е) - £-1*J. (14)

Вычислим внешний дифференциал от уравнения (14):

д

dt(dE) = d(£-1*dH) - d(£-1*(a*E)) - d(£-1*J). (15)

Рассмотрим уравнение (15). Для перехода к уравнению второго порядка необходимо получить дифференциал в явном виде. Распишем выражение d(£-1*(a*E)) с целью выделения дифференциала dE:

d(£-1*(a*E)) = d(£-1*a*E) = d(£-1aE). (16)

Выражение под дифференциалом в (16) является произведением тензоров, и по свойствам дифференцирования тензорных величин имеем [10]

d(AB) = (dA) B + A (dB). (17)

По свойствам тензора (17) перепишем (15) в виде

д

— (dE) = d(£-1*dH) - d(£-1a)E - £-1adE - d(£-1*J). (18)

Полагая, что каждая компонента тензора не зависит от пространственных координат, подставим (1) в (18):

д2 д

d(£-1*dH) + — (H) + £-1а—(^*H) - d(£-1*J) = 0. (19)

д 2 д

Вывод 2. Вычисление внешнего дифференциала от уравнения (2) сразу, без умножения на е-1 и оператор Hodge *, приводит к неправильным результатам, поскольку получаемый член уравнения d(dH) обращается в нуль по свойствам дифференциальных форм. Переход к уравнению второго порядка для магнитного поля накладывает ограничения на компоненты тензора электропроводности. Компоненты тензора могут быть либо константами, либо функциями только времени.

4. Система уравнений Максвелла в анизотропных средах

Запишем систему уравнений Максвелла, используя аппарат дифференциальных форм с учетом материальных соотношений (3), (4):

др*Н

dt

-dE;

де*Е dt

-- dH - а*Е - J;

dD = p; dB = 0.

(20)

(21)

(22) (23)

Здесь плотность заряда р — дифференциальная форма порядка 3 (3-форма).

Таким образом, из уравнений (20)-(23) видно, что все константы, входящие в систему уравнений Максвелла, могут быть тензорами ранга 2. В частности, электрическая проводимость а является тензором ранга два.

Для выполнения необходимых условий непрерывности компонент поля (13) на межфрагментарной границе подобластей с различными значениями тензорной электропроводности требуется выполнение слабой дивергенции

div

д

ш (еЕ) + аЕ

0

или в терминах дифференциальных форм

d

д

- (е*Е) + а*Е

0.

(24)

Докажем, что для уравнений Максвелла первого порядка условие (24) выполняется автоматически.

Лемма 1. Доказать, что для любых ф Е Н*(П), ф — 0-форма,

' д (е*Е)

d

dt

+ а*Е

Л фdQ = 0.

По свойствам дифференциальных форм преобразуем интеграл к виду

d

д(е*Е) dt

+ а*Е

а

~Id

п

Л фdQ = ( д(е*Е)

V dt

д (е*Е)

+ а*Е Л ф

+ а*Е

dQ.

Л dфdQ-

(25)

п

п

п

Второе уравнение Максвелла в терминах дифференциальных форм выглядит следующим образом:

d

д(е*Е) ,

dt

+ а*Е

Л WdQ = J (|-1*Б) Л dWdQ+

п , п , (26) + d((/j,-1 *Б) Л W)dQ+ J Л dфdП,

п п

где W — тестовая функция из пространства H(rot; П).

Согласно диаграмме де Рама dф является 1-формой и принадлежит пространству H(rot;Q), тогда, применяя (26) к (25), получим

" д (е*Е)

d

дt

+ а*Е

Л фdП = J (|-1*Б) Л d^)dQ + J J Л dфdП+ пп

+ J d(^-1*B Л dф)dn + J d^

- а*Е) Л ф^П.

(27)

п п

По свойствам дифференциальных форм ¿(¿ф) = 0 и

^ 3 Л ¿фдП = J ¿3 Л ф^ + у ¿(3 Л ф)^. п п п

Условием разрешимости уравнений Максвелла является условие ¿3 = 0. С учетом этого (27) принимает вид

" д (£*Е)

d

дt

+ а*Е

Л фdП = J d(/i-1*B Л dф)dП+ п

+ / d ( ( -- а*Е + J ) Л ф ) dQ.

(28)

д (е*Е)

Из второго уравнения Максвелла следует, что ——--+ а*Е- J = d(|-1*Б). Тогда (28)

можно записать в виде " д (е*Е)

дt

d

at

+ а*Е

Л фdП = d [(|-1*Б) Л dф] dQ - d [d(/i-1*B) Л ф] dQ. (29)

Перепишем (29) " д (е*Е)

d

дt

+ а*Е

Л фdП= d(|-1*Б) Л dфdП+ (|-1*Б) Л d(dф)dП-

- d ^(|-1*Б)] Л фdП - d(|-1*Б) Л dфdП = 0.

(30)

пп Лемма доказана. □

Вывод 3. Для уравнений Максвелла первого порядка показано автоматическое выполнение условий непрерывности компонент поля на межфрагментарной границе в анизотропной области. Система уравнений Максвелла первого порядка не накладывает каких-либо ограничений на компоненты тензоров параметров среды £, ^, а. С учетом вышеизложенного при моделировании электромагнитных полей в анизотропных средах целесообразнее решать уравнения Максвелла как систему уравнений первого порядка.

п

п

п

п

п

п

п

п

п

5. Математическая модель

Источники возмущения электромагнитного поля могут быть как искусственными (поле проводников с током), так и естественными (естественное электромагнитное поле Земли). Моделирование электромагнитного поля в данной работе основано на декомпозиции поля на первичное и вторичное поля [11]. Первичное поле — падающая ТЕ (ТМ)-волна, вторичное поле — отраженная волна, полученная от взаимодействия падающей волны с проводящим полупространством.

Имеет место следующая декомпозиция поля:

E = Ep + Es, B = Bp + Bs.

Подставив (31) в систему уравнений Максвелла, получим

д Bs / dBp

- rot Es - —- + rot EM ,

dt \ dt dEs s _ ( дEp

е— = rot /i-1Bs - oEs - е— - rot /i-1Bp + oEP

H0(div; П) = {v е H(div; П) : rot v = 0} , H00(div; П) = {v е H0(div; П) : v ■ n|r = 0}

(31)

dt ^ V dt

Обозначим

dBp

Jm = ~gf + rot Ep; (32)

д E

Je = е— - rot ^-1BP + oEP, (33)

e dt K '

где Jm и Je — распределенные источники поля.

Первичное поле в виде ТЕ (ТМ)-волны имеет следующий вид:

TE : E =(Ex(x,y,z),Ey(x,y,z), 0), H = (0, 0,HZ(x,y,z)),

TM : E = (0, 0, Ez(x,y,z)), H = (Hx(x,y, z), Hy(x,y,z), 0).

6. Вариационная формулировка

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

H(П) = {v : v е [L2(П)]3} ,

H(rot; П) = |v е ^2(П)]3 : rot v е ^2(П)]3} , H0(rot; П) = {v е H(rot; П) : div v = 0} , H(div; П) = {v е ^2(П)]3 : divv е L2^)} , H00(rot; П) = {v е H0(rot; П) : v x n|r = 0} ,

где n — вектор внешней нормали к границе Г области Q. Здесь L2(Q) — пространство интегрируемых с квадратом функций. Векторные элементы могут быть разделены на два класса элементов: конформных в пространстве Н(div; Q) и конформных в Н(rot; Q) [2, 3]. Для пространства Н(rot; Q) вводят скалярное произведение

(u, v)H(rot;n) = J u ' vdQ + J rot u ■ rot vdQ.

пп

Скалярное произведение для пространства Н(div; Q) вводят следующим образом:

(u, v)H(div;n) = J u ■ vdQ + J div u ■ div vdQ.

пп

Для введенных пространств выполняются условия включения:

Vp e L2(Q) grad p e Н(rot; Q),

VE e Н(rot; Q) rot E e Н(div; Q). (34)

Тогда можно сформулировать следующую задачу:

Для заданных источников поля Jm и Je найти Es e H0(rot; Q) и Bs e Н(div; Q) такие, что:

д Bs

- rot Es - Jm в Q; (35)

е— = p-1 rot Bs - aEs - Je в Q; (36)

дИ

д Es

"at

Es x n|ri = Er; (37)

Bs ■ n|r2 = Br, (38)

где Г = Г1 U Г2.

Введем скалярное произведение

(u, v)= uvdQ.

Выпишем вариационную постановку для задачи (35)-(38). Для этого в соответствии с введенным скалярным произведением уравнение (35) умножим на функцию из пространства ^(div, Q) и уравнение (36) — на функцию из пространства Н°(гс^; Q). Тогда получим следующую вариационную постановку.

Найти Es e Н°(го^ Q) и Bs e Н(div; Q) такие, что для VW e Н°(rot;Q), VF e Н00(div; Q)

' д Bs \

д ,F) = -(rot Es, F) - (Jm, F),

. дt дEs \

е-Qj-, W j = (p-1 rot Bs, W) - (aEs, W) - (Je, W). Используя векторное тождество

div(a x b) = b ■ (rot a) - a ■ (rot b),

п

получим

J rot y-1B ■ Wdn = J rot W ■ ^-1BdH + J div(W х y-1B)dП. (39)

n n n

Применяя формулу Гаусса — Остроградского

У divu dn = У (u ■ n)dS,

П дП

получим

-1 -1

j div(W х ^-1B)dH = J (W х y-1B) ■ ndr. (40)

n Г

Из определения тестовой функции W следует, что

J (W х y-1B) ■ ndr = 0. (41)

Г1

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

Найти Es G H0(rot; П) и Bs G H(div; П) такие, что для VW G H(rot;Q), VF G H0(div; П)

Bs ■ Fdn = J rot Es ■ F -J Jm ■ F в П; (42)

nnn

дJ £E ■ Wdn = J y-1 rot W ■ Bdn -J aE ■ WdH -J Je ■ WdH+

n n (43)

. — 1тэ s\

+ J (W х y-1Bs) ■ ndr в П.

Г2

7. Дискретизация вариационной постановки

Пусть в области П выполнено разбиение на множество согласованных параллелепипедов. Определим дискретное пространство Hh(rot; П) С Hr(rot;n), введя векторные базисные edge-функции на одном элементе. Дискретное пространство Hh(div; П) С HT(div; П) определим, введя векторные базисные face-функции на одном элементе.

Используя дискретные пространства Hh(rot; П) и H^div;^, сформулируем дискретный аналог вариационной постановки (42), (43).

Найти Es G Hh(rot; П) и Bs G Hh(div; П) такие, что для VW G H^rot;^, VF G Hh(div; П)

д J Bs ■ FdП = J rot Es ■ F -J Jm ■ F на П,

n n n

дJ £E ■ WdП = J y-1 rot W ■ BdП -J aE ■ WdП -J Je ■ WdП+ n n n n

+ У (W х y-1Bs) ■ ndr на П.

Г2

Разложим искомые функции Е и В по базисам пространств Ин(тоЬ; П) и Ин(&у; П):

Е = £ егтг, В = £ Ъ3(1)¥3.

3

Поскольку Jm, Je — распределенные источники, их можно также представить в виде разложения по базисам:

Je = ^ зер(г)~пр, Jm = ^ ]тк(г)¥к. р к

Здесь Нр — edge-функции, определенные на ребрах параллелепипеда; Гк — Гаев-функции, определенные на гранях параллелепипеда.

В результате получим систему обыкновенных дифференциальных уравнений относительно коэффициентов ег(г) и Ъз (г):

СЗЪ(г) = -Ке() - Рт, Сде(г) = КТЪз (г) - 8ег(г) - + Z.

Используя для аппроксимации по времени соотношения

де еп+1 - еп

дг т

дЪ Ъп+1 - Ъп-2

дг т

получим систему линейных алгебраических уравнений:

ОЪп+2 = ОЪп-2 - тКеп - тР]тп+ 2,

Сеп+1 + т8еп+1 = Сеп + тКТЪп+2 - тЦзеп+2 + тZn+1,

(44)

где элементы матриц О, С, К, 8, Р, Я и вектора правой части Z определяются соотношениями

[0]г,з = I /1-1¥г • ¥3вП; (45)

п

[С\3 = ! ^ • N3вП; (46)

п

[8]г,з = ! aN, • N3вП; (47)

п

[Р]г,з = ! • ¥звП; (48)

п

[Я]гз = I N • N3вП; (49)

п

[К]г,з = ! /л-1¥г • N3вП; (50)

п

[Z ]i

(Ni х ц-1В) ■ ndr.

(51)

Г2

Здесь (45)-(49) — матрицы массы, а матрица (50) является матрицей произведения функций edge- и face-базисов. В результате вычислений получим конкретный вид локальных матриц в изотропной среде для базиса в пространстве H(div, П)

P _ lxlylz

2 1 0 0 0 0

1 2 0 0 0 0

0 0 2 1 0 0

0 0 1 2 0 0

0 0 0 0 2 1

0 0 0 0 1 2

и базиса в пространстве H(rot, П)

Q

Qxx

0 0

0 о

Qyy 0

0Q

Qxx = Qyy = Q

Матрица (50) имеет вид

lxlylz

36

4 2 2 1 2 4 12 2 14 2 12 2 4

(K)

T

-1

2lxlz lxlz -2lxly -l l xy 0 0

2l xl z -l l xz -l l xy -2lxly 0 0

lxlz 2lxlz 2lxly lxly 0 0

-l l xz -2lxlz lxly 2lxly 0 0

2ly lz -l l yz 0 0 2lxly lxly

-l l yz -2lylz 0 0 2lxly -l l xy

2lylz lylz 0 0 lxly 2lxly

lylz 2l l yz 0 0 -l l xy -2lxly

0 0 2l l yz lylz -2lxly -l l xy

0 0 -2lylz -l l yz -l l xy -2lxly

0 0 lylz 2lylz 2lxly lxly

0 0 -l l yz 2ly lz lxly 2lxly

Определим конкретный вид локальной матрицы для тензорного коэффициента а

с XX аху $2

S

&yxS2 Oyy S1 Oyz S2 &zx S3 & zy St Ozz. Si

S1

lxlylz

36

zy 2

'zz S1

4221 2412 2142 1224

6

1х1у

24

1х1у

24

2 12 1 2 12 1 12 12 12 12

2 2 11 1122 2211 1122

8. Тензорный коэффициент электрической проводимости в анизотропных средах

В изотропных средах электрическая проводимость а является скаляром или диагональной матрицей 3 х 3 с равными значениями на диагонали. В случае трансверсально-изотропной среды электрическая проводимость изменяется вдоль одной из осей (рис. 1), а остается диагональной матрицей, но значение ах отличается от ах и ау в Л2 раз:

а

Ох

О О

О

ау

О

О О

Ох

ахх О О

О ауу О

О О ахх

1 < Л2 < 2.

В случае наклонно-изотропной среды электрическая проводимость изменяется слоями, расположенными под углом к поверхности земли (рис. 2). Введем новую систему координат х',у',2', начало координат которой совпадает с началом координат исходной системы х,у,2 и оси повернуты на угол —в. В новой системе координат среда является трансверсально-изотропной и тензор а является диагональным. Выпишем выражение для а в исходной системе координат:

а(х,у,2) = 3а(х' ,у' ,2),

(52)

где 3 — матрица перехода от системы х', у',2' в систему х,у, 2.

Получим уравнения осей координатной системы х', у' ,2' в системе х,у ,2. Спроектируем точку А, лежащую на оси О'х' и имеющую координаты (1,0,0) в системе х',у',2' на

£21, еО, цО, <т1=10е-6

£21, еО, /лО, о1=10е-6

О х

Рис. 1. Трансверсально-изотропная среда. Рис. 2. Наклонно-изотропная среда.

оси системы x,y,z. Получим точку A c координатами (cos 9, 0, — sin 9) в системе x,y,z. Аналогично спроектируем точку B', лежащую на оси O'y' и имеющую координаты (0,1,0) в системе x',y', z' на оси системы x, y, z. Получим точку B c координатами (0, cos 9, — sin 9) в системе x, y, z. Тогда, применив формулу для уравнения прямой, проходящей через две точки к парам точек: 0(0, 0,0) и A, 0(0, 0, 0) и B, получим уравнения для осей O'x' и O'y' в координатах системы x,y,z соответственно:

O'x' :

x

cos 9

sin 9

O y :

y

cos 9

sin 9

(53)

Через точки O, A, B и оси O'x', O'y' проходит плоскость, описываемая уравнением

sin 9x + sin 9y + cos 9z = 0. (54)

Поскольку ось O' z' является перпендикуляром к плоскости, описываемой (54), не труд-

но выписать ее уравнение:

O z : —

x

x

sin 9

sin 9 cos 9

(55)

Тогда из (53), (55) следует, что матрица J принимает вид

J

cos 9 0 — sin 9 0 cos 9 — sin 9 — sin 9 — sin 9 cos 9

(56)

Подставив (56) в (52), получим выражение для электрической проводимости наклонно-

изотропной среды:

a

ax cos 9 0

0

—az sin 9

ay cos 9 —az sin 9 —ax sin 9 —ay sin 9 az cos 9

9. Вычислительный эксперимент

Тест 1. Моделируется вторичное поле над полупространством в области П. Первичное поле имеет форму ТЕ-волны вида

Ex = cos nx sin ny cos(wt — z),

Ey = — sin nx cos ny cos(wt — z),

П

Hz = — 2— cos nx cos ny(sin(wt — z) + sin(z)). w

Тогда токи Jm, Je, определенные соотношениями (32), (33), принимают вид

— sin nx cos ny sin(wt — z) Jm =| — cos nx sin ny sin(wt — z)

2n(1 + cos nx cos ny cos(wt — z)

^ — cos nx sin ny[—ew sin(wt — z) —— (sin(wt — z)+ sin z) + a cos(wt — z)] ^

w^

Je = 2n2

sin nx cos ny [—ew sin(wt — z)--(sin(wt — z) + sin z) + a cos(wt — z)]

w^ 0

z

z

z

Q2, (72, £0, f.10

Z //

У

z / £2i, ai, so, /ло /

Рис. 3. Расчетная область.

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

1

(рис. 3), £ = £0, Ц = Ц0, а,

воздух

10

а

-6

Ом

м

аземля — тензор вида

ахх 0 0 аУУ (

0 0 а

где ахх = аУУ = 0.1, azz = 0.5ахх, или

а

ахх cos в 0

—ахх sin в — аУУ sin в

—azz sin в

аУУ cos в —azz sin в

zz zz

azz cos в

где в = 0, 30,60, 90°.

В табл. 2 и 3 представлены значения компонент электрического поля и магнитной индукции в точке (0.4, 0.45, 0.51) на 1000-м шаге по времени. На рис. 4 показаны силовые линии в горизонтальном сечении z = 0.51.

Тест 2. Моделируется вторичное поле над полупространством в области П. Первичное поле имеет вид ТМ-волны:

п

Нх = — sin пх cos ny(sin(wt — z) + sin(z)), w п

НУ =--cos пх sin ny(sin(wt — z) + sin(z)),

w

Ez = sin пх sin пу cos(wt — z). Тогда токи Jm,Je, определенные соотношениями (32), (33), принимают

п(1 + ц) sin пх cos пу cos(wt — z) —п(1 — ц) cos пх sin пу cos(wt — z) 0

Jm

п

Je

--cos пх sin ^(cos z — cos(wt — z))

w

--sin пх cos пу (cos z — cos(wt — z))

w

2п2

, sin пх sin пу(а cos(wt — z) — £w sin(wt — z)--(sin(wt — z) + sin z) ,

w

0

Таблица 2. Зависимость значений компонент электрического поля в точке (0.4,0.45,0.51) от свойств среды

Однородная среда Диагональный тензор

Ех -0.0580066 -0.0580066

Еу 0.0282856 0.0282856

Ех 2.17965в-5 2.18716в-5

Плотный тензор(0 = 0°) Плотный тензор(0 = 30°)

Ех -0.0580066 -0.0596304

Еу 0.0282856 0.0291886

Ех 2.18716в-5 -0.000915992

Плотный тензор(0 = 60°) Плотный тензор(0 = 90°)

Ех -0.0639195 -0.0700522

Еу 0.0313926 0.0344734

Ег -0.00175316 -0.0022833

Таблица 3. Зависимость значений компонент магнитной индукции в точке (0.4,0.45,0.51) от свойств среды

Однородная среда Диагональный тензор

Вх 6.54175Б-11 6.5415Б-11

Ву 1.3426Б-10 1.34257Б-10

Вг 4.38117Б-11 4.38117Б-11

Плотный тензор(0 = 0°) Плотный тензор(0 = 30°)

Вх 6.5415Б-11 8.71919Б-11

Ву 1.34257Б-10 1.55305Б-10

Вг 4.38117Б-11 4.2346Б-11

Плотный тензор(0 = 60°) Плотный тензор(0 = 90°)

Вх 1.08122Б-10 1.23219Б-10

Ву 1.76853Б-10 1.93851Б-10

Вг 3.82352Б-11 3.22215Б-11

а б

Рис. 4. Горизонтальные компоненты электрического поля (Ех,Еу) (а) и магнитной индукции (Вх,Ву) (б) в сечении г = 0.51.

Таблица 4. Зависимость значений компонент электрического поля в точке (0.4,0.45,0.51) от свойств среды

Однородная среда Диагональный тензор

ЕХ -0.00120045 -0.00144901

ЕУ -0.000589735 -0.000711668

Ех 0.139202 0.158243

Плотный тензор(0 = 0°) Плотный тензор(0 = 30°)

ЕХ -0.00144901 0.00532095

ЕУ -0.000711668 0.00610844

Ех 0.158243 0.163221

Плотный тензор(0 = 60°) Плотный тензор(0 = 90°)

ЕХ 0.0122319 0.0178147

ЕУ 0.0131437 0.0188998

Ех 0.177546 0.197017

Таблица 5. Зависимость значений компонент магнитной индукции в точке (0.4,0.45,0.51) от свойств среды

Однородная среда Диагональный тензор

ВХ -6.93195Б-11 -7.83089Б-11

ВУ 1,41312Б-10 1.59662Б-10

Вг 6.11882Б-16 7.4863Б-16

Плотный тензор(0 = 0°) Плотный тензор(0 = 30°)

ВХ -7.83089Б-11 -6.37479Б-11

ВУ 1.59662Б-10 1.47578Б-10

Вг 7.4863Б-16 -3.37745Б-12

Плотный тензор(0 = 60°) Плотный тензор(0 = 90°)

ВХ -5.4183Б-11 -5.149144Б-11

ВУ 1.45044Б-10 1.5185Б-10

Вг -6.87381Б-12 -9.73255Б-12

а б

Рис. 5. Горизонтальные компоненты электрического поля (ЕХ,ЕУ) (а) и магнитной индукции (ВХ,ВУ) (б) в сечении г = 0.51.

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

1

£ — £q , ^ — ^q, ^воздух

10

-6

, аземля тензор вида

м

ахх 00

а = 0 ayy 0 ,

0 0 azz

где ахх — а.

yy

0.1, аг

0.5ах

<т

х, или

ахх cos 9 0

—azz sin 9

ayy cos 9 —azz sin 9

—ахх sin 9 — ayy sin 9

'zz 1 ?zz '

azz cos 9

где 9 — 0, 30, 60, 90 В табл. 4 и 5 представлены значения компонент электрического поля и магнитной индукции в точке (0.4, 0.45, 0.51) на 1000-м шаге по времени. На рис. 5 показаны силовые линии в горизонтальном сечении z — 0.51.

Список литературы

[1] Hiptmair R. Multilevel Preconditioning for Mixed Problems in Three Dimensions: Dissertation zur Erlangung des Doktorgrades der Mathematisch. Naturwissenschaftlichen Fakultät at der Universität at Augsburg, Augsburg, 1996. 145 p.

[2] Nedelec J.C. Mixed finite elements in R3 // Numer. Mathematik. 1980. Vol. 35, N 3. P. 315341.

[3] Nedelec J.C. A new family of mixed finite elements in R3 // Numer. Mathematik. 1986. Vol. 50. P. 57-81.

[4] Raviart A., Thomas J.M. A mixed finite element method for 2-d order elliptic problems. Mathematical aspects of finite elements // Proc. of the Conf. Held in Rome, 10-12, Dec. 1975 / A. Dold, B. Eckman (Eds). Heidelberg, N.Y.: Spinger-Verlag, Berlin, 1977. P. 292-315 (Lecture Notes in Mathematics. Vol. 606).

[5] Bossavit A. Differential Geometry for the Students of Numerical Methods in Electromagnetism. Electricite de France, Etudes et Recherchers, Aug., 1991.

[6] White D.A. Discrete Time Vector Finite Element Methods for Solving Maxwell's Equations on 3D Unstructured Grids: PhD Thesis, Lawrence Livermore National Laboratory, Sept., 1997. 245 p.

[7] Gradinarü V., Hiptmair R. Multigrid for discrete deferential forms on sparse grids // Computing. 2003. Vol. 71. P. 17-42.

[8] Абрамов А.А. Введение в тензорный анализ и риманову геометрию. М.: Физматлит, 2004.

[9] Hiptmair R. Discrete Hodge-operators: an algebraic perspective // Progress in Electromagnetics Research, PIER 32. 2001. P. 247-269.

[10] Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. М.: Наука, 1973.

[11] Aruliah D.A. Fast Solvers for Time-Harmonic Maxwell's Equations in 3D: PhD Thesis, The University of British Columbia, Aug., 2001. 160 p.

0

Поступила в редакцию 1 марта 2006 г.