Расчёт эффективных электрофизических характеристик в многомасштабной изотропной среде Текст научной статьи по специальности «Математика»

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

Том 17, № 6, 2012

Расчёт эффективных электрофизических характеристик в многомасштабной изотропной среде

Е.П. Курочкинл1, О.Н. Соболева2 1 Институт теплофизики СО РАН, Новосибирск, Россия 2Институт вычислительной математики и математической геофизики СО РАН,

Новосибирск, Россия e-mail: kurochkina@itp.nsc.ru, olga@nmsf.sscc.ru

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

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

Введение

В геофизических задачах крупные неоднородные включения (пласты, пропластки) учитываются в математической модели непосредственно с помощью граничных условий (см., например, [1-3]). Пространственное распределение мелкомасштабных неоднородностей редко известно точно и часто описывается случайными полями. Поэтому задачи для сред с вариациями физических параметров на всех масштабах требуют значительных вычислительных затрат. Традиционный подход к решению задач, включающих малые масштабы, состоит в поиске более простых моделей, требующих меньших вычислительных затрат, решение которых для физических величин, например, напряжённости электрического поля, плотности тока, было бы близко в среднем к решению первоначальной полной задачи. Построение таких моделей, правильно описывающих поведение решения в крупномасштабном пределе, в литературе известно как гомогенизация, огрубление сеток, подсеточное моделирование, расчёт эффективных коэффициентов [4-8]. Эти подходы наиболее развиты в теории стационарной фильтрации [6-10].

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

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

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

1. Постановка задачи

Согласно [14] уравнения Максвелла для монохроматических полей Е (х, ¿) = И,е х х (Е (х) е-^), Н (х,г) = Яе (Н (х) е-^) имеют вид

гоШ (х) = (—Ъше(х) + а (х)) Е (х) + Е,

го1Е = гш^Н, (1)

где Е и Н — векторы напряжённости электрического и магнитного полей; е (х) — диэлектрическая проницаемость; ^ — магнитная проницаемость; а (х) — электропроводность; ш — циклическая частота; х — вектор пространственных координат. Магнитная проницаемость ^ равна магнитной проницаемости вакуума. В неограниченной области выполняются условия излучения от источника Е, т.е. на бесконечности решение системы уравнений (1) затухает. При а(х)/ (ше(х)) ^ 1, когда токи проводимости значительно преобладают над токами смещения, диэлектрическая проницаемость среды е (х) слабо влияет на характеристики поля; амплитуда и фаза поля зависят в основном от электропроводности среды а (х). В этом случае задачу рассматривают в квазистационарном приближении. При больших сопротивлениях среды на высоких частотах появляется зависимость измеряемого сигнала от диэлектрической проницаемости. Циклическая частота, электропроводность и диэлектрическая проницаемость удовлетворяют неравенству а(х)/ (ше(х)) < 1.

Для моделирования полей а(х), е(х) используется подход, подробно описанный в работе [15]. Пусть поле электропроводности а(х) известно. Это означает, что выполнено его измерение на некотором масштабе /0 в каждой точке х. Чтобы перейти к более грубой сетке масштабов, недостаточно сгладить а (х)г по масштабу I, I > 10, так как сглаженное поле не будет правильно отражать физический процесс, описываемый уравнениями (1) на интервале масштабов (1,Ь), где Ь — максимальный масштаб неоднородности среды. Это объясняется тем, что флуктуации проводимости на интервале масштабов (/0,/) коррелируют с флуктуациями напряжённости электрического поля Е и эти корреляции могут быть достаточно большими.

Для построения модели среды, как и в работе [16], рассматривается безразмерное поле ф, равное отношению полей, полученных сглаживанием проводимости а (х)г по двум различным близким к (/0) масштабам 1,1'. Обозначим через а(х) сглаженное по масштабу I поле а1о (х). Тогда ф(х, I, I') = а(х)^/а(х), I' < I. Случайное поле ф меняется плавно по сравнению с полями а(х)^/, а(х)г. Раскладывая поле ф в ряд относительно I' — I и оставляя только члены первого порядка малости при I' ^ I, получим уравнение

д 1п а(х) .

^Пг = (2)

где <^(х,/') = (дф(х, I', 1'у)/ду) |у=1. Фактически мелкомасштабные флуктуации поля (р могут наблюдаться только в некотором конечном диапазоне масштабов 10 < I < Ь. Решение уравнения (2) имеет вид

а1о(х) = ао ехр | — / ^(хЛ)^ | , (3)

1о

где о0 — константа. Согласно теореме о суммах независимых случайных полей [18], если дисперсия <^(х, /) в данной точке конечна, то при больших значениях ¿//0 интеграл в (3) стремится к полю с нормальным распределением вероятностей. Если дисперсия поля <^(х, /) бесконечна и существует невырожденное (не сосредоточенное в одной точке) предельное распределение суммы случайных величин, то это распределение является устойчивым. В данной работе предполагается, что поле <^(х, /) имеет нормальное распределение и изотропную однородную корреляционную функцию

< р(х,/) р(у,/') > - < р(х,/) >< р(у,/') >=

= Фдд(|х - у | ,/,/')£ (1п/ - 1п /') (4)

(угловые скобки означают статистическое усреднение). Из формулы (4) следует, что флуктуации поля ^ в разных масштабах не коррелируют. Это обычное предположение для скейлинговых моделей соответствует тому факту, что статистическая зависимость становится незначительной в случае, если масштабы флуктуаций параметров различны по величине [16]. Если же среда масштабноинвариантная, то для любого положительного значения К выполняется условие

Фдд(|х - у | , /, /') = Фдд(К |х - у| , К/, К/').

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

£10(х) = ео ехр х(хЛ)^

V 1о

Предполагается, что функция х(х,/) имеет нормальное распределение вероятностей и дельтокоррелированна по логарифму от масштаба /. Если диэлектрическая проницаемость для любого / удовлетворяет условию < еДх) >= е0 и масштабноинвариантна, то

ФХХ = 2 < X > . (6)

Корреляционная функция между полями проводимости и диэлектрической проницаемостью определяется корреляцией между полями х(х,/) и <^(х,/):

Фдх(х,у,/,/') = (р(х, /)х(у, /')> - (р(х, /)>(х(у, /')> = Фдх(|х - у| , / , ВД1п/ - 1п/'). (7)

При не масштабноинвариантной среде величины Фдд, Фо%, Фдх зависят от масштаба /.

2. Подсеточная модель

Функции проводимости и диэлектрической проницаемости о(х) = о(х)г0, е(х) = е(х)г0 разделим на две компоненты относительно масштаба /. Крупномасштабные (надсеточ-ные) компоненты о (х, /), е (х, /) получены статистическим усреднением по всем /1) и х(х,/1) для /0 < /1 < /, / - /0 = где мало. Мелкомасштабные (подсеточные)

компоненты равны а'(х) = а(х) — а(х, /), е'(х) = е(х) — е(х, /):

е(х, /) = е0 ехр

ехр

. ¿0

е'(х) = е(х,/)

ехр

— I I1

ехр

—IК*«^

-1

(е'(х)) = 0,

а(х, /) = а0 ехр

— / ^(хЛ) ТГ

ехр

— / ^(хЛ) ТГ

а'(х) = а(х,/)

ехр

— / ^(хЛ) 17

. ¿0

ехр

—/ ^(хл) ТГ

10

1

(а'(х)) = 0.

Из формул (8) следует, что с точностью до членов второго порядка малости

е(х, /) ~ а(х, /) ~

1 1

1/

1 — (х) т- + офХх (/) т

/ 2 1/

1 — Мт + (От

/2

(x), а1(х).

(9)

Крупномасштабные компоненты напряжённости электрического и магнитного полей Е (х,/), Н (х,/) получаются как усреднённые решения системы уравнений (1), в которых крупномасштабные компоненты а(х,/), е (х, /) фиксированы, а мелкомасштабные а'(х), е'(х) — случайные поля. Подсеточные компоненты электрического и магнитного полей равны Е' (х) = Е (х) — Е (х, /), Н' (х) = Н (х) — Н (х,/). Подставим выражения для Н (х), Е (х), а(х), е (х) в систему уравнений (1) и усредним по мелкомасштабным компонентам:

гоШ (х, /) = (—¿ше (х, /) + а (х, /)) Е (х, /) + ((—¿ше' + а') Е') + Е,

:ю)

го1Е (х, /) = ^¿шН (х, /).

Подсеточный член ((—гше' + а') Е') в системе (10) не известен и не может не учитываться без предварительной оценки. Несмотря на то что мелкомасштабные компоненты е', а' малы, корреляция с подсеточной напряжённостью электрического поля может быть значительной. Оценка этого члена определяет подсеточную модель. Подсеточный член оценивается с помощью теории возмущений. Вычитая систему (10) из системы (1) и оставляя только члены первого порядка малости, получим подсеточные уравнения

гоШ' = (—¿ше (х, /) + а (х, /)) Е' + (—гше' (х) + а' (х)) Е (х, /)

го1Е' = ^¿шН'.

ь

Переменная состояния Е (х, /) в правой части (11) считается известной. Решение системы уравнений (11) равно [19]

1 [ 1

Е' =—гшЛ -е}кг (—гше' (х') + а' (х')) Еа (х',/) Лх'+ 4п ] г

1 д д С 1

+ 4п (—гше (х,0 + а (х,/)) дХадХв / 1 ^ (х'> + а'(х'» Е ^" (12)

где г = |х — х'|, к2 = ш^. (ше (х, /) + га (х, /)). Для определённости выбрано то значение корня, при котором И,е к > 0, 1т к > 0. Из (12) следует, что подсеточный член в (10) равен

<(—гше' (х) + а' (х)) Е' (х)> = = -1 гш^J1 егкг <(—гше' (х) + а' (х)) (—гше' (х') + а' (х'))> Еа (х', /) Лх'+

, / (—гше' (х) + а' (х)) д д [1 . , , , , \ , ,

+ ( . , .-- / ,чч „ _егкг ( —гше' (х') + а' (х')) ) Ее (х', /) Лх'. (13)

\4п (—гше (х, /) + а (х,/)) джа джв У г у у ' у '7 в ^ у

Для полей, в которых небольшое изменение масштаба влечёт значительные изменения самого поля (это характерно для сильно меняющихся физических параметров, описываемых мультикапликативными каскадами), можно считать, что а(х, /), е (х, /), Е (х, /) и их производные меняются медленнее, чем а', е', Е' и их производные. Поэтому а(х, /), е (х, /), Ее (х', /) можно выносить за знак интеграла. Переходя к сферическим координатам, интегрируя по частям и оставляя только первые члены малости Л///, при условии а(х)/ (ше(х)) < 1 получим оценку

<(—гше' (х) + а' (х)) Е'а (х)> «

те

« — 1 (2^ш2е(х, /) — гш^а (х, гегкгФхх (г) Лг^гше(х, /)Еа (х, /) +

о

те

+ 2 (2^ш2е(х, /) — гш^а (х, /)) ^ гегкгФх<7 (г) ^а(х, /)Еа (х, /) +

о

те

Г Л/

+гш^а(х, /) / гегкгФст<7 (г) Лг^-а (х, /) Еа (х, /) +

о

+1ФХХ у гше(х, /)Е (х, /) + — ^о)) уа(х, /)Ег (х, /). (14)

Если ш^Ь2 |(гше (х,/) + а (х, /))| ^ 1, то интегралы в (14) малы [20]. Поскольку максимальный масштаб неоднородностей много меньше длины волны, это неравенство не является слишком ограничительным. Можно записать

1 Л/

<—гше' (х) Е'а (х)> + <а' (х) Е'а (х)> « — ^ (0) (—гше(х, /)Ев (х, /)) у —

— ((0) — 1фХХ (0^ !/а(х, /)Еа (х, /). (15)

Подставим (15) в (10):

гоШ (х, /) = —гше10 ехр

аю

е10

— / х(х,/1)*

Е (х,/) + аю ехр

— / ^(х,11) 17

го1Е (х, /) = гш^Н (х, /) ,

1

ФХХ 1/ ТГ 7

1 +

¥—и) I'

е0,

1— ( 2фХ° (0) — 1фХХ (о^ у

1 +

ао.

Е (х, /)

16)

Из (16) с точностью до членов второго порядка по 1/// следует, что новые коэффициенты а^0 и ею удовлетворяют следующим равенствам:

е10 = е0 +

ф*х , Л 1/ 1Т — (х)ео7,

1

1

а,о = ао + — -ФХ0 + тгФХХ + 7^0° — (?) а„^.

3

3

2

1/

/

Устремляя в этих равенствах 1/ к нулю, получим уравнения

11п ео1 1

т^Х* — (х)

6

11п аог 11п /

11п / 211 —^Х0 + ^ФХХ + ^0° — М.

:17)

В масштабноинвариантной среде решение данной системы уравнений зависит от масштаба / степенным образом и эффективные уравнения имеют вид

гоШ (х, /) = — ег (х) Е (х, /)

(х)-фХх/6

/ \ <0>) + 2ФХ^-3ФХХ-1

аг (х) Е (х,/) ,

го1Е (х, /) = гш^Н (х, /)

18)

3. Численное моделирование

Для проверки приведённых выше формул численно решается задача (1). Используются

следующие безразмерные переменные: х = х/—0, а = а/а0, Н = Н/Н0, ЕЕ

-оао

кНО

Е,

ше0 ао

. В расчётах к = 5, к1 = 4л/2. Таким образом,

к1 = —0^а0^ш, к = к^а — гке, к задача решается при а0 = 1, е0 = 1. В безразмерном виде уравнения имеют вид

гоШ = (а — гка) к1ЕЕ, гс^Е = гкьН.

19)

Область интегрирования разбивается на три подобласти. В области 0 < а < 0.1, 1.5 < а < 1.6 помещается поглощающий слой [21]. Для этого производится замена

ь

ь

переменных (Жг) =

1 + гп(Жг), которая обеспечивает экспоненциальное затухание волны в поглощающем слое. Функция п(Жг) для наших данных выбиралась в виде

0.1 - Жг'

п(Жг)

3.5 0,

3.5

0.1

1.6 — Жг

0.1

0 < Жг < 0.1, 1.1 < Жг < 1.5, 1.5 < Жг < 1.6.

Таким образом, переменные £г отличаются от переменных Жг только в поглощающем слое. В областях 0.1 < Жг < 0.3, 1.3 < Жг < 1.5 коэффициенты Х и у равны 1. В точке (0, 0, 0.2) находится источник ^ = 0, = 0, = 0.5ехр — д2(Ж3 — 0.2)2, д = 60. В области 0.3 < Жг < 1.3 проводимость и диэлектрическая проницаемость моделируются мультипликативными каскадами. Интегралы (3), (5) аппроксимируются конечно-разностными формулами, в которых удобно перейти к логарифму по основанию 2:

а (х)

¿0

е (х)1

ехр

ехр

— 1п 2

^ (х, тг) дьт

1°§2 ¿0

— 1п 2

X (х т)

1°§2 ¿0 1, I = 2Т, Ат

- 12 ^(х,Т4)аТ

2 4=-8 ,

- 12 х(х,Т4)аТ

2 4=-8 .

(20)

Здесь (а(х)1о) = 1, (е(Х)1о) = 1, I = 2Т, Ат — шаг по т. В расчётах Ат выбирается равным 1. В показателе степени (20) стоит сумма статистически независимых слагаемых, поскольку предполагается статистическая независимость для полей с различными масштабами. Для расчётов используются два слагаемых с г = —5, г = —4. Остальные слагаемые полагаются равными нулю. Количество слагаемых выбиралось так, чтобы масштаб самых крупных вариаций проводимости и диэлектрической проницаемости позволил заменить приближённо вероятностные средние величины осреднёнными по пространству, а самых мелких — так, чтобы разностная задача хорошо аппроксимировала уравнения (19). При каждом г поля строятся по формулам

р(х, тк)

Ш2

(1(Х,тк) +

Фхх / \

х(Х,тк) = у (гС1(Х,тк) + v/Г-72C2(X,тfc)) +

фхх

2 V ' ,24 ' "V 2 '

где поля ^1, С2 — независимые, гауссовы, с единичной дисперсией, нулевым средним и корреляционной функцией

«1(Х, тг)С1(у, т,-))с = «2(Ж, тг)С2(у, т,-))с = ехр [— (х — у)2 /22т*¿г,] .

Коэффициент ФдХ в этом случае равен г^/ФХхФ0°. Поля С2 моделируются с помощью метода, изложенного в [22]. При расчётах г полагалось равным единице. Коэффициенты Фо должны выбираться из экспериментальных данных. В нашем расчёте они равны 0.4, а < ^ >=< х >= 0.2. На рис. 1 приведено поле проводимости для выбранных

2

о о

Рис. 1. Поле проводимости для двух масштабов г = —5, —4 вычисленное по формуле (20) в среднем сечении хз = 1/2

а

0.4 0.6 0.8 1 1.2 х3 ' 0.4 0.6 0.8 1 1.2

Рис. 2. Результаты расчёта компонент реальной и мнимой частей напряжённостей магнитного поля по оси х2 (а) и электрического поля по оси х\ (б): 1 — для усреднённых коэффициентов (при а = 1, £ = 1); 2 — для эффективных коэффициентов; 3 — усреднённый по 48 реализациям результат, полученный с помощью прямого численного моделирования

масштабов в среднем сечении x3 = 0.5. Для решения уравнений (19) использовались метод, основанный на конечно-разностной схеме Yee, и метод декомпозиции, предложенный в [23]. Последний позволяет разбить задачу на восемь подзадач, каждая из которых решалась независимо на параллельной машине. Затем собиралось общее решение. Задача решалась на сетке 412 х 412 х 412 с постоянным шагом по всем переменным 1/256. Результаты расчётов приведены на рис. 2. В соответствии с процедурой вывода подсеточных формул для проверки необходимо много раз решить полную задачу и выполнить вероятностное усреднение по мелкомасштабным вариациям (кривая 3). В итоге получится решение, которое можно сопоставить с решением эффективных уравнений (18). Для сравнения приводится также решение уравнений со средними значениями коэффициентов проводимости и диэлектрической проводимости (кривая 1). Усреднение по ансамблю Гиббса требует многократного решения полной задачи. В работе применяется более экономный вариант проверки. Решение при каждом значении x3 усредняется по плоскостям (xi, x2), а затем проводится доусреднение по ансамблю. В расчётах использовались 48 реализаций.

Результаты численного моделирования показали, что длина волны для усреднённых полей увеличивается на 5 % по сравнению с длиной, полученной для средних коэффициентов (кривые 1, 3, см. рис. 2). Эффективные коэффициенты (кривая 2) дают ошибку по длине волны только 0.5 %. Ошибка по амплитуде для эффективного решения меньше в два раза по сравнению с ошибкой, которую дают усреднённые коэффициенты.

Заключение

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

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

[1] Mikhailenko B.G., Soboleva O.N. Mathematical modeling of seismomagnetic effects arising in the seismic wave motion in Earthes's constant magnetic field // Appl. Math. Lett. 1997. Vol. 10, No. 3. P. 47-55.

[2] Млстрюков А.Ф., МихАйлЕнко Б.Г. Численное решение уравнений Максвелла в анизотропных средах на основе спектрального преобразования Лагерра // Геология и геофизика. 2008. № 8. С. 819-829.

[3] Эпов М.И., Шурина Э.П., Нечаев О.В. Прямое трёхмерное моделирование векторного поля для задач электромагнитного каротажа // Там же. 2007. № 9. С. 989-995.

[4] Yukalov V.I., Glüzman S. Self-semilar bootstrap of divergent series // Phys. Rev. E. 1997. Vol. 55. P. 6552-6570.

[5] Yukalov V.I. Self-semilar approximations for strongly interacting systems // Phys. A. 1990. Vol. 167. P. 833-860.

[6] Gluzman S., Sornette D. Self-similar approximants of the permeability in heterogeneous porous media from moment equation expansions // Transport in Porous Media. 2008. Vol. 71. P. 75-97.

[7] Germano M., Moin P. Piomelly U., Cabot W.H. A dynamic subgrid scale eddy viscosity model // Phys. Fluids. A. 1991. Vol. 3, No. 7. P. 1760-1765.

[8] Germano M., Sagat P. Large Eddy Simmulation for Incompressible Flow. Berlin, Heidelberg: Springer, 1998.

[9] Sahimi M. Flow phenomena in rocks: From continuum models, to fractals, percolation, cellular automata, and simulated annealing // Rev. of Modern Phys. 1993. Vol. 65. P. 1393-1534.

[10] Dagan G. Higher-order correction of effective permeability of heterogeneous isotropic formations of lognormal conductivity distribution // Transport in Porous Media. 1993. Vol. 12. P. 279-290.

[11] Biswal B, Manwart C., Hilfer R. Three-dimensional local porosity analysis of porous media // Phys. A. 1998. Vol. 255. P. 221-241.

[12] Бобров Н.Ю., Лювчич В.А., Крылов С.С. Масштабная зависимость кажущегося сопротивления и фрактальная структура железистых кварцитов // Изв. РАН. Физика Земли. 2002. № 12. C. 14-21.

[13] Bekele A., Hudnall H.W., Daigle J.J. et al. Scale dependent variability of soil electrical conductivity by indirect measures of soil properties //J. of Terramech. 2005. Vol. 42. P. 339-351.

[14] Ландау Л.Д., Лифшиц Е.М. Электродинамика сплошных сред. М.: Наука, 1982.

[15] Кузьмин Г.А., Соболева О.Н. Подсеточное моделирование фильтрации в пористых автомодельных средах // ПМТФ. 2002. Т. 43, № 4. C. 115-126.

[16] Kolmogorov A.N. A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number //J. Fluid Mech. 1962. Vol. 13. P. 82-85.

[17] Монин А.С., Яглом А.М. Статистическая гидромеханика. Т. 2. М.: Наука, 1975.

[18] Гнеденко Б.В., Колмогоров А.Н. Предельные распределения для сумм независимых случайных величин. Л.: Гостехиздат, 1954.

[19] Глинер Э.Б., Кошляков Н.С., Смирнов М.М. Основные дифференциальные уравнения матаматической физики. М.: Физматгиз, 1962.

[20] Кравцов Ю.А., Рытов С.М., Татарский В.И. Введение в статистическую радиофизику. М.: Наука, 1978.

[21] Chew W., Jin J., Michielssen E., Song J. Fast and Efficient Algorithms in Computational Electromagnetics. Artech House, 2001.

[22] Ogorodnikov V.A., Prigarin S.M. Numerical Modeling of Random Processes and Fields: Algorithms and Applications. Utrecht, Netherlands, 1996.

[23] Davydycheva S., Drushkin V., Habashy T. An efficient finite-difference scheme for electromalnetic logging in 3D anisotropic inhomogeneous media // Geophysics. 2003. Vol. 68. P. 1525-1536.

Поступила в 'редакцию 21 августа 2012 г.