Многосеточный метод решения прямой задачи индукционного каротажа Текст научной статьи по специальности «Математика»

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

Том 13, Специальный выпуск 4, 2008

Многосеточный метод решения прямой задачи индукционного каротажа

Н. А. Луговая

Омский государственный университет им. Ф.М. Достоевского, Россия

e-mail: lugovaya@omsu.ru

The problem of induction logging is reduced to a boundary value problem in a domain containing both a well bore and the invaded zone. We use an integro-differential boundary condition for the case of a horizontally layered containing medium. We present a multigrid iterative method for solution of a two-dimensional axisymmetric problem and provide the results of numerical experiments for test problems.

Введение

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

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

Рассмотрим изотропную среду с электромагнитными параметрами: диэлектрической проницаемостью е, электрической проводимостью о и постоянной магнитной проницаемостью, равной проницаемости вакуума, ^ = Электромагнитное поле изменяется во времени по закону е~гшЬ с круговой частотой и.

Рассмотрим области П С П2 С Я3, Область П с границей Г1 и внешней нормалью п1 является частью области скважины. Она содержит источники электромагнитного поля и имеет постоянные электромагнитные параметры. Данные источники в однородном пространстве Я3 с параметрами среды П1 возбуждают первичное электромагнитное поле Е0, Н0. Область П2 с границей Г2 и внешней нормалыо п2 имеет скважину и зону проникновения. Электромагнитные параметры внешней к П2 области Д3\П2 — функции только одной декартовой координаты (например, в горизонтально-слоистой среде).

Предполагается, что комплексная проводимость о' = о — 1ие является кусочно-непрерывной и ограниченной функцией, В точках непрерывности о' электромагнитное поле описывается уравнениями Максвелла, На поверхностях разрыва о' выполняются

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

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

Будем рассматривать электрическую составляющую Е электромагнитного поля, а в области ^ уравнения будем записывать относительно аномального поля Е — Е0,

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

Рассмотрим гильбертово пространство

Н(г<Л, П) = {и| и е (£2(П))3, г<Ли е (Ь2(П))3},

||и|| = J {|и|2 + |го!и|2} СП.

п

Вариационная постановка задачи имеет вид: найти Е е Н(го1;, П2) такое, что для любого

V е Н(г<^, П2)

У {гоt(E — Е0)ШУ — к2(Е — Е0)У} СП + J {г<ЛЕг<^V — к2ЕУ} СП =

П1 П2\П1

= — у [г<ЛЕ X У]п2СГ ^ у [п^Е0 х У^СГ; (1)

Г2 Г1

■[п2 х то1Е] - \п2 х [{и ц[п2 х --[та2 х го^пйС -

Г2

2'

-^(гг2го^)а1уСя}СГ =0, к2 = ги/мт'. (2)

Интегродифференциальное граничное условие (2) следует из обобщения формул Стрэт-тона—Чу для случая горизонтально-слоистой вмещающей среды [2]. При этом Он — матрица фундаментальных решений магнитного типа для задачи в пространстве К3 с параметрами вмещающей среды Д3\П2, Интегрирование в уравнении (2) следует понимать в смысле главного значения Коши, поскольку диагональные элементы ОНР(М, М0)

при совпадении точек М и М0 имеют особенность типа —^——.

КММо

Рассмотрим случай вертикальной скважины, когда свойства среды не зависят от азимутальной координаты <р (0 < ^ < 2п) цилиндрической системы координат {г, г}, (г > 0, —то < г < оо), ось Ог которой совпадает с осью скважины, В качестве источника выберем набор соосных со скважиной катушек, возбуждающих только азимутальную компоненту электрического поля Е^, Тогда задача решается в пространстве П1 С П2 С К+ х К с искомой функцией и(г, г) = Е^(г, г) и известным первичным полем и0 (г, г) = Е° (г, г) в гильбертовом пространстве:

П(то^, П) = {и I и е ь2(П), ^ е ь2(П), е ь2(П)

дг г дг

и =

1 д(ги) т дг

т йт ¿г.

При этом решение удовлетворяет условиям:

1)и{0,г) — и0{0, г) = 0 У г е {—ж, ж);

2) и ~ о(1/К), —--гки ~ о(1/К) при К = л/г2 + г2 —> сю,

дт

Будем также предполагать, что область ^ совпадает с областью скважины, а контур Г2 является вертикальной прямой. Тогда вариационная постановка (2), (2) имеет вид У у е Н{гоtlp, П2)

ду 1 д д(ту)

д_

дг'

1д_

т дт

П2\П1

ди ду 1 д{ти) д{ту) , 2

'7Г7Г + —\ \ ~ к гиу дг дг т дт дт

1 д , ч П, 2 „ д

дт йтйг =

д{ти) дт

Г2

2 дт

(ти) +

Г2

к гиС + — (ги) —--г

дт дт

дС дидС

дг дг

уйг — J

Г1

¿г = 0.

йтйг +

' д(ги°) дт

у ¿г;

(3)

(4)

Фундаментальное решение С{М, М0) данной задачи в пространетве Я+ х Я с электромагнитными параметрами вмещающей среды Я+ х К \ представляется потенциалом кольцевого пояса вертикально ориентированных магнитных диполей [3]. При

совпадении точек М и М0 ядро С(М, М0) имеет логарифмическую, а ядро --^——

дг

сильную особенности, поэтому интегрирование последнего слагаемого следует понимать в смысле главного значения Коши,

3. Дискретизация

Решение задачи (3), (4) будем искать в прямоугольнике Б = {{т,г) | 0 < т < т2, г1 < г < г2}, граница которого выбирается согласно условиям:

1) {г,т2) е Г2 У г е {—ж, ж);

2) и{т, г) ж 0 Ут > 0, У г е {—ж, г-\) и {г2, ж).

Область покроем прямоугольной сеткой с шагом к = {кг ,кг), Определим множества узлов сетки

Б1 = {{г,зЖтг,г3) е Б}, Г = {{г,з)1{тг,г3) е Г2}

и дискретное подпространство Нн{^2) = врап{ е Н{го1^, П2), (г,з) е Б1}, а также подмножества

О1! = {{г,3) е Б1 I вирр Ыз П = 0}, В* = {{г,з) е Б11 вирр игз П П2\^1 = 0}. Представив решение в виде линейной комбинации и кг е Нн{Б)

и{т, г) — и0{т, г) ж ^ {и кг — ик^и^т, г), {т,г) е Пь (к,1)ео'1

и(

ь(г,г) ~ ^ иыУы{г,г), (г, г) Е

выписываем (3) для всех {итп1(т,п) Е Он} и приходим к системе уравнений

дукг дьтп , 1 д(гукг) д(гьтп)

^ (ик1 - и°к1)

ик1 /

дг дг г дг дг дькг дьтп . 1 д(гькь) д(гУтп)

- к гУкхУтп

дг дг

Г д(ги) дг

г дг дг д(ги0)

к гик1итп

¿гё,г =

дг

(5)

Г2

Г1

В работе в качестве базисных функций итп использовались билинейные функции, равные 1 в (т, п)-узле и 0 в остальных. Интегралы по контуру Г2 аппроксимировались с четвертым порядком точности.

Для дискретизации (4) воспользуемся методом коллокаций. Решение на границе представляем в виде

и

Г2

£

кеГ

®к Vк,

д(ги)

дг

Г2

Рк ук,

(6)

кеГ

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

ак = \uk-i ~ Юик + ик+х} , (Зк =

д(ги)

- 10

д(ги)

д(ги)

дг ) к-1 V дг )к V дг /к+1 Выписывая (4) во всех точках коллокации М) Е Г", приходим к системе уравнений

(7)

+

д(ги) дг

+ V Ли = 0,

(8)

где

д(ги) дг

и и — сеточные функции, определенные в узлах из Г^. Матрица Л опреде-

ляСТСя формулами (7), , сс дИаГ„„алЫ1Ие алсмсптъ. равПИ 12, а эломоПТИ „ад . „од

8

диагональю равны —. Элементы матриц ¿Т н V:

8

к

д

дг

■С(М,М1)ик(зм) Мм,

Г2

V)

1к

д

дУк

к2гС(М,М1)ик(зм) -^—С(М,М1)-^(зм)

дг

дг

¿£л

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

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

д(ги)

Выражая из системы (8) значения функции —-— и подставляя их в правую часть

дг

системы (5), приходим к системе уравнений с неизвестными значениями функции и на ОН: 1и + Zu = /.Здесь 1 — оператор левой части равенства (5), а оператор Z действует только на граничные точки контура Г2,

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

Многосеточный метод представляет собой итерационный процесс, на каждой итерации которого приближенное решение задачи ищется за счет сглаживания невязки и рекурсивного спуска па более редкую сетку [1]. В связи с этим через Ф обозначим пространство функций иН, определенных на сетке ОН, Для простоты изложения предположим, что Нг = 2-псг и Нг = 2-тсг.

При квазистационарной постановке задачи (Б,е к2 = 0) в качестве оператора, эффективно гасящего негладкую компоненту невязки, использовался метод Зейделя, Для восстановления значений сеточных функций на более частую сетку использовался оператор линейной интерполяции Рк : Ф

,1, -^2/г, -т-И

м • гь _ ф ; такой, что

иН

К ,2Н

г]

Р ки

и

2к

г] '

0.5(и2-

1]

и

2 Н )

г+1] ),

2 Н

0-5(и2]-1 + и2]+1),

0-25(и2—1]-1 + и2+1]-1 + и2-1]+1 + и2+1]+1)

если «^даны, если г нечетно и ] четно, если г четно и нечетно, если г и нечетны.

—к

—2/1

А при переходе на более редкую сетку использовался оператор Цк : Ф" —> сопряженный к РН и определяемый равенством

и

2 Н

г]

Яки] =0,2иНк1 + 0.5(и] + ^ иНк1) ( |к - г|< 1, 11 — ]|< 1 ). (9)

(к-г)(1-]) = 0

Систему, построенную на сетке ОН, запишем в виде

Ькик = / Н.

(10)

Одна итерация многосеточного метода состоит в следующем. Сделав в ~ 4 итераций методом Зейделя, имеем приближение иН, Из равепства ЬНиН — /к = г к определяем невязку г Н на сетке Ок. Решив задачу вида Ькик = г Н, мы найдем поправку ик к приближенному решению ик и определим точное решение задачи (10) как ик = ик — ик.

Поскольку невязка в данном случае является гладкой функцией, то с небольшой погрешностью поправку можно найти, решая аналогичную задачу на более редкой сетке, Для этого строим оператор Ь2Н на сетке 02Н. Значения интегралов, составляющих оператор </, пересчитываем по рекурсивным формулам, а оператор Z перестраивается согласно формуле (9), Преобразовав систему к виду

и решив ее приближенно, определяем приближенное значение и)2Н. Взяв новое приближенное решение йн = йн — Рсглаживаем певязку в итерациями метода Зейделя с противоположным порядком перебора узлов.

Для решения задачи (11) па сетке 02Н используется та же самая процедура с пулевым начальным приближением и вспомогательной сеткой 04Н. Таким образом, рекурсивный спуск продолжается до тех пор, пока на очередной вспомогательной сетке число узлов не станет совсем незначительным, И на этой редкой сетке система легко решается прямым методом.

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

5. Численный эксперимент

Результаты работы многосеточного метода покажем на следующей модельной задаче, В качестве источника выбирается вертикально ориентированный магнитный диполь на оси скважины при г = —0.1 м. Радиус скважины полагается равным 0.08 м, Рассматривается цилиндрически слоистая среда с электромагнитными параметрами:

ь2Нт2Н = ф V

(Н)

0.8

о 0.6

о

а

а

и 0.4

о

0

о

О2 0.4 0.6

0.8

г, м

Погрешность , -—

\ЕЧ>{г г)\

% на отрезке г € [0.1, 1.1м], г = 0.03 м при вычислениях на сетках 26 х 29 и 27 х 210

В прямоугольнике О = [0,1.28] х [—5, 5] задача решается на двух сетках размерностью 27 х 29 и 28 х 210, На обеих сетках сделано четыре итерации многосеточного метода, и по окончании невязка алгебраической системы составляла порядка 10-7,

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

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

[1] Дьяконов Е.Г. Минимизация вычислительной работы. М.: Наука, 1989. 272 с.

[2] Дмитриев В.И., Захаров Е.В. Интегральные уравнения в краевых задачах электродинамики. М.: МГУ, 1987.

[3] Табаровский Л.А. Применение метода интегральных уравнений в задачах геоэлектрики. Новосибирск: Наука, 1975.

Поступила в редакцию 20 февраля 2008 г.