CC BY
9
УДК 519.632.4+550.832.7
DOI: 10.183 03/2618-981X-2018-4-162-170
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ СИГНАЛОВ ЭЛЕКТРОМАГНИТНОГО ЗОНДА С ТОРОИДАЛЬНЫМИ КАТУШКАМИ В ДВУМЕРНЫХ ИЗОТРОПНЫХ МОДЕЛЯХ ГЕОЛОГИЧЕСКИХ СРЕД
Ирина Владимировна Суродина
Институт вычислительной математики и математической геофизики СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Лаврентьева, 6, кандидат физико-математических наук, старший научный сотрудник; Институт нефтегазовой геологии и геофизики им. А. А. Тро-фимука СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Коптюга, 3, старший научный сотрудник, тел. (913)903-87-41, e-mail: sur@ommfao1.sscc.ru
В работе рассматривается математическое моделирование электромагнитных откликов нового каротажного прибора в осесимметричных изотропных средах, выполненное с помощью конечно-разностного метода. Описываются постановка задачи, методы решения. Приводятся численные эксперименты, позволяющие оценить величину влияния корпуса прибора на измеряемые сигналы в зависимости от длины корпуса. Обсуждаются результаты моделирования.
Ключевые слова: уравнение Гельмгольца, каротаж, математическое моделирование, конечные разности, система линейных алгебраических уравнений, итерационные методы, прямые методы.
NUMERICAL SIMULATION OF ELECTROMAGNETIC TOOL SIGNALS WITH TOROIDAL COILS IN TWO-DIMENSIONAL ISOTROPIC GEOLOGICAL MEDIA MODELS
Irina V. Surodina
Институт вычислительной математики и математической геофизики СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Лаврентьева, 6, Ph. D., Senior Researcher; Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 3, Prospect Аkademik Koptyug St., Novosibirsk, 630090, Russia, Senior Researcher, phone: (913)903-87-41, e-mail: sur@ommfao1.sscc.ru
In this paper we propose the mathematical simulation of electromagnetic responses of the new logging tool in axially symmetric isotropic media, to be carried out with the help of a finite difference method. To solve the resulting system of equations, we use iterative methods with implementation on the CPU, GPU and a direct method (Intel MKL PARDISO). A comparative analysis is given. The results of numerical experiments are presented and discussed. The simulation results were used to justify the optimal configuration when developing a new device.
Key words: Helmgolz equation, logging, simulation, finite differences, system of linear algebraic equations, iterative methods, direct methods.
Введение
В последнее десятилетие произошли большие инженерные достижения в каротажном приборостроении. Появилась возможность создавать приборы, основанные на нетрадиционных для каротажа источниках возбуждения элек-
тромагнитного поля. В ИНГиГ и НПО «Луч» создан прибор нового поколения -устройство для регистрации характеристик электромагнитного поля с использованием тороидальных катушек [1-7]. В процессе создания прибора применялось математическое моделирование. Это самый быстрый и дешевый способ воспроизвести работу зонда в моделях, близких к реальности. Целью настоящего исследования было правильно сформулировать задачу для осесимметричных моделей, найти метод решения и оценить влияние величины корпуса прибора на измеряемые характеристики. Результаты данного исследования были учтены в [3, 4] при создании оптимальной конфигурации зонда [5].
Постановка задачи
Рассмотрим случай изотропной проводящей немагнитной среды. Для амплитуд электрического и магнитного поля справедливы уравнения Максвелла:
гоШ = аЕ ШЁ = т\\Н - у
(1)
_п
где а - электропроводность, = ]ы0 = 471-10 Гн/м - магнитная проницаемость, со - циклическая частота, у - источник стороннего тока (плотность тока,
Е - электрическое поле, Н - магнитное поле). Зависимость тока от времени гармоническая. Будем рассматривать только осесимметричные модели сред. Введем цилиндрическую систему координат {г, ф, 2}. Учитывая тот факт, что в осесимметричной среде плотность тока имеет только одну ненулевую компоненту - тангенциальную ]0 = {0, ,0}, = -г0Жг ~го) ( КР-
магнитный момент, 8 - дельта функция Дирака), система уравнений (1) упрощается и принимает вид:
дИ ф
1Иф+дИф
г дг
дЕ дЕ„
<зЕ„
(2)
д2 дг
т^И,
ф
Iц
ф
Выразим компоненты электрического поля через компоненту магнитного поля и подставим в последнее уравнение (2):
<
<
д_
дг
1 дИ,
Ф
чст дг у
д_
дг
д ( гИ, ) уог дг у
(3)
Из симметрии модели и источника следует условие на оси скважины:
Иф 1г=0 = 0 •
(4)
В соответствии с условием излучения на бесконечности магнитное поле
затухает с удалением от источника
И.
ф
г, г ^<х>
->0. Это позволяет приближен-
но поставить нулевые граничные условия для компоненты И, на большом расстоянии от источника:
Иф 1г =Я = 0, Иф \г=±2 = 0 •
(5)
Вместе с условием (4) получаем для магнитного поля разностную краевую задачу Дирихле (3)-(5) для уравнения Гельмгольца (3).
Для построения конечно-разностной схемы [8] источник расположим в точку (= 0, г0). В силу осевой симметрии рассмотрим полуплоскость [0, Я] х [-2, 2], где введем прямоугольную неравномерную координатную сетку
= {(г-,X I = а.^Ыг,] = -Ы2^^Ы2}, N = Я, = 2 •
(6)
Рассмотрим на сетке (6) линейное конечномерное пространство Н0 сеточных функций, равных нулю на границе. Скалярное произведение в Н0 определяется следующим образом:
Ыг -1 Ы2 -1
(7)
г=0 )=1-Ы2
(г)
где щ' - средний шаг на неравномерной сетке по переменной г.
= + 2, =Ц~ Ц-\, г = 1,... Средний шаг по переменной ъ
определяется аналогично. и, V - сеточные функции. На пространстве Н0 определим разностный оператор А
АУ = -(ЬУ?)г -(«с(гУ)г)я-дУ
Л
где V е И0, с, q - сеточные функции, с = (г - кк-'Г / 2) 1, д = /юц, Ь, а - сеточные функции сопротивления среды (р = -1), которые в узлах двумерной сета
ки (¡, у) определены так:
к кц к кц
Ьц = Р(г + у,2 ц- у), аи = Р(г - у,2 ц + у).
Для разностных производных на неравномерной сетке (6) приняты обозначения [9]:
(У)г,и = <?и -= (Ум] -Уф/(Пг,и = (У1+и -Уф//£>.
Разностные производные по переменной z определяются аналогично. После дискретизации получим разностное уравнение
AV = F, (8)
где /¡ц = - цф. Упорядочим двумерные векторы V и F по столбцам, трансформируя их в одномерные массивы. Тогда (8) будет представлять систему линейных алгебраических уравнений с пятидиагональной блочной матрицей А. Полученная матрица является комплексной, не эрмитовой, но симметричной
в пространстве Н0 со скалярным произведением (7). В пространстве, где скалярное произведение определяется без участия шагов сетки, данная матрица будет не симметричной. Для уменьшения количества операций при решении уравнения (8) перейдем в пространство с обычным скалярным произведением
Ыг-1 Ыг-1
(= I I иг]Уг] (9)
/=0 ]=1-Ы2
и симметризуем матрицу А с помощью диагонального преобразования подобия [10] А = Б~1/2А01/2, 0 = ^(4,...,^), е Яп.
В итоге получим систему линейных алгебраических уравнений
АХ = Е\ (10)
где Х = Б"1/2У, ¥ = В1/2¥
Численные эксперименты
Для решения полученной системы с симметричной матрицей применялись итерационные методы и прямой метод. Первые расчеты были проведены с по-
165
мощью метода эрмитова разложения [11], в дальнейшем он заменен методом сопряженных A-ортогональных сопряженных невязок (Conjúgate A-Ortogonal Conjúgate Residual (COCR)) [12] как более подходящий для данного класса задач и требующий меньшее количество итераций. Итерационные методы были реализованы с предобуславливателем из метода последовательной симметричной верхней релаксации (SSOR). Программы написаны автором. Реализация прямого метода осуществлена с помощью программы PARDISO из библиотеки INTEL MKL. Итерационный процесс завершался при получении решения заданной точности (относительно точного решения). Как правило, относительное падение невязки составляло 8 порядков. И прямой метод, и итерационные показали хорошие результаты по точности получаемых решений. Быстродействие будет обсуждено ниже в данной статье. Сравнение разных решателей подтверждает правильность работы программ.
Первые эксперименты были проведены для моделей, имеющих точное решение. Это цилиндрические вертикально-слоистые и горизонтально-слоистые изотропные среды. По результатам этих экспериментов были отобраны сетки, погрешность решения на которых не превышает 2 % для достаточно большого диапазона сопротивлений в моделях (от 0,02 до 500 Омм). Расчеты проведены для следующих частот: 5, 25, 50, 100 КГц. Сравнивались расчетные величины в приемных катушках от 0,2 до 1,1 м.
Для оценки влияния корпуса зонда использована следующая модель.
Прибор: радиус 0,51 м, сопротивление корпуса 5,7*10 "8 Ом м.
Среда: скважина радиуса 0,108 м имеет сопротивление 1 Омм, пласт характеризуется сопротивлением 10 Омм.
Длина металлического корпуса варьировалась от 1 до 100 м. На рис. 1 приведены расчеты для реальной и мнимой части компоненты в точке приемной катушки - 0,2 м и точное решение для бесконечного корпуса для прибора с частотой 5 КГц.
Рассмотрим рис. 1. По оси абсцисс отложена полудлина металлического корпуса в метрах. По оси ординат приведены значения вещественной и мнимой
частей компоненты , рассчитанной для корпуса зонда с конечными длинами
(от 1 до 100 м), и точные значения Яф (бесконечный корпус).
Как видно из рисунка, для частоты 5 КГц выход на асимптоту (точное решение) для вещественной и мнимой частей компоненты происходит только
при длине корпуса, большей 80 м (полудлине 40 м). Это позволяет сделать следующий вывод. Для моделирования тороидального зонда (учитывая реальные размеры корпуса) с низкими частотами необходимо двумерное моделирование. Асимптотикой в данном случае пользоваться нельзя, слишком велико влияние размера корпуса. Аналогичная картина получается и для зонда с более высокой частотой - 100 КГц, расчеты для которого в приемной катушке 0,2 м приведены на рис. 2 (для приемных катушек 0,4-1,1 м результаты аналогичны).
Рис. 1. Зависимость от длины корпуса. Зонд с частотой 5 КГц:
сплошные линии - реальная часть И, с положительными значениями и мнимая
часть И, с отрицательными значениями; соответствующие асимптоты (прерывистые линии) для бесконечного корпуса
Рис. 2. Зависимость от длины корпуса. Зонд с частотой 100 Кгц:
сплошные линии - реальная часть И, с положительными значениями и мнимая
часть И, с отрицательными значениями; соответствующие асимптоты (прерывистые линии) для бесконечного корпуса
Рассмотрим рис. 2. Здесь также по оси абцисс отложена полудлина корпуса. Отметим 2 факта. Первый: для частоты 100 Кгц выход на асимптоту наблюдается значительно раньше (18-20 м), чем для низкочастотного зонда. Второй: абсолютные значения компонент в 6 раз больше, чем для зонда с частотой 5 КГц.
Из приведенных выше расчетов следует, что для моделирования сигналов тороидальных зондов (с частотами от 5 до 100 КГц) необходимо двумерное моделирование (в осесимметричных средах), поскольку достаточно велика разница между расчетными сигналами в радиально-слоистой среде с бесконечным корпусом и возможным реальным прибором.
Вычисление компоненты Ez и обсуждение результатов
Одна из задач данной работы заключалась в исследовании поведения сигналов в зависимости от длины корпуса. Это было сделано для компоненты .
На практике наряду с измерениями тангенциальной компоненты магнитного поля измеряется еще и вертикальная компонента электрического поля. Эту компоненту также необходимо уметь точно рассчитывать. На первый взгляд, это не должно составить особого труда - компонента определяется явно из второго уравнения системы (2). Однако при расчетах пришлось столкнуться с некоторыми проблемами при вычислении производной. Это связано с очень высоким контрастом сопротивлений между корпусом и скважиной (до 9 порядков), поведением компоненты Ez - мнимая часть имеет большой градиент и меняет знак в измеряемом диапазоне (в приемных катушках от 0,2 до 1,1 м).
При вычитании компонент происходила потеря значащих цифр и в результате получалась низкая точность расчетной компоненты, не удовлетворяющая заданному диапазону в 2 %. Эта проблема была решена с помощью использования свойства непрерывности этой компоненты и очень мелкого шага в узлах,
участвующих в вычислении односторонней производной - 7.8* 10_6.
Использование такого мелкого шага катастрофически увеличило время вычислений с использованием итерационного метода - от нескольких секунд
(для сеток, когда вычислялась только компонента ) до нескольких часов.
Тогда расчеты были перенесены на графический процессор. Написаны специальные программы на основе полностью параллельных алгоритмов, аналогично [13]. На GPU был реализован метод COCR с предобусловливателем, предложенным нами в [14]. Время вычислений сократилось в 80 раз, но все еще оставалось достаточно большим (несколько минут).
Использование мелкого шага для расчетов, где реализован прямой метод (PARDISO), лишь незначительно (порядка 3 %) увеличило время вычислений, которое составило несколько секунд. Таким образом, для вычислений электрической и магнитной компонент предпочтительнее использовать прямой метод, реализованный в библиотеке INTEL MKL.
Заключение
Создана программа прямого двумерного моделирования для зонда с тороидальными катушками. Проведены тестовые расчеты на правильность работы программы. Показана зависимость показаний зонда от длины корпуса. Результаты данных расчетов были приняты во внимание при изготовлении прибора.
Благодарности
Автор выражает признательность коллегам по лаборатории скважинной геофизики Института нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН за содержательное обсуждение постановки задачи и благодарит В. Н. Глинских за постоянные консультации, М. Н. Никитенко за предоставленную программу расчета тороидального зонда в слоистой среде, М. И. Эпова за полезные советы и замечания.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Способ измерения удельной электропроводности и электрической макроанизотропии горных пород: пат. РФ № 2525149 / М. И. Эпов, В. Н. Глинских, М .Н. Никитенко ; опубл. 10.08.2014 ; бюл. № 22.
2. Устройство для измерения удельной электропроводности и электрической макроанизотропии горных пород: пат. РФ № 2528276 / М. И. Эпов, В. Н. Еремин, А. К. Манштейн и др. ; опубл. 10.09.2014 ; бюл. № 25.
3. Устройство для регистрации характеристик электромагнитного поля с использованием тороидальных катушек: пат. РФ № 2578774 / М. И. Эпов, В. Н. Еремин, А. Н. Петров и др. ; заявл. от 14.01.2015 ; опубл. 27.03.2016 ; бюл. № 9.
4. Устройство для генерации электромагнитного поля тороидальной катушкой в геологической среде: пат. № 2579177 / М. И. Эпов, В. Н. Еремин, А. Н. Петров и др. ; заявл. от 14.01.2015 ; опубл. 10.04.2016 ; бюл. № 10.
5. Электромагнитный зонд для каротажа в нефтегазовых скважинах: пат. РФ № 2583867 / М. И. Эпов, В. Н. Еремин, А. Н. Петров, В. Н. Глинских ; опубл. 10.05.2016 ; бюл. № 13.
6. Математическое и физическое моделирование сигналов электромагнитного зонда для изучения макроанизотропии осадочных отложений [Электронный ресурс] / М. И. Эпов, В. Н. Глинских, В. Н. Еремин и др. // Геомодель 2017: 19-я конференция по вопросам геологоразведки и разработки месторождений нефти и газа - Геленджик, 11-14 сентября 2017 г.: тезисы докладов. (43809). - URL: http://www.earthdoc.org/publication/publicationdetails/ ?publication=90886 (дата обращения - 01.03.2018).
7. Electromagnetic tool for high-resolution logging: theoretical and experimental studies [Electronic source] / M. I. Epov, V. N. Glinskikh, V. N. Eremin et al. // SPE Russian Petroleum Technology Conference. Moscow, October 16-18, 2017. - P. SPE-187904-MS. - URL: https://www.onepetro.org/conference-paper/SPE-187904-MS (дата обращения - 01.03.2018).
8. Мартаков С. В., Эпов М. И. Прямые двумерные задачи электромагнитного каротажа // Геология и геофизика. - 1999. - Т. 40, № 2. - С. 249-254.
9. Самарский А. А. Теория разностных схем. - М. : Наука, 1989. - 616 с.
10. Кузнецов Ю. И., Агапитова Н. С. Математические основы моделирования на ЭВМ. - Ю.-Сахалинск : Из-во ЮСИЭПИ, 2003. - 214 c.
11. Воеводин В. В., Кузнецов Ю. А. Матрицы и вычисления. - М. : Наука,1984. - 318 c.
12. Sogabe T., Zhang S. A COCR method for solving complex symmetric linear systems // Journal of Computational and Applied Mathematics. - 2007. - Vol. 199. - P. 297-303.
13. Surodina I. The GPU Solvers for High-Frequency Induction Logging // Numerical Analysis and Its Applications, NAA 2016, Lecture Notes in Computer Science. - Vol. 10187. -P. 640-647. - URL: https://link.springer.com/content/pdf/10.1007 %2F978-3-319-57099-0.pdf
14. Labutun I. B., Surodina I. V. Algorithm for Sparse Approximate Inverse Preconditioners in Conjugate Gradient Method // Reliable Computing (Interval Computations) Journal. - 2013. -URL: http://interval.louisiana.edu/reliable-computing-journal/tables-of-contents.html#Volume_19.
REFERENCES
1. Sposob izmerenija udel'noj jelektroprovodnosti i jelektricheskoj makroanizotropii gornyh porod: pat. RF № 2525149 / M. I. Jepov, V. N. Glinskih, M. N. Nikitenko ; opubl. 10.08.2014 ; bjul. № 22.
2. Ustrojstvo dlja izmerenija udel'noj jelektroprovodnosti i jelektricheskoj makroanizotropii gornyh porod: pat. RF № 2528276 / M. I. Jepov, V. N. Eremin, A. K. Manshtejn i dr. ; opubl.
10.09.2014 ; bjul. № 25.
3. Ustrojstvo dlja registracii harakteristik jelektromagnitnogo polja s ispol'zovaniem toroidal'nyh katushek: pat. RF № 2578774 / M. I. Jepov, V. N. Eremin, A. N. Petrov i dr. ; zajavl. ot
14.01.2015 ; opubl. 27.03.2016 ; bjul. № 9.
4. Ustrojstvo dlja generacii jelektromagnitnogo polja toroidal'noj katushkoj v geologicheskoj srede: pat. № 2579177 / M. I. Jepov, V. N. Eremin, A. N. Petrov i dr. ; zajavl. ot 14.01.2015 ; opubl. 10.04.2016 ; bjul. № 10.
5. Jelektromagnitnyj zond dlja karotazha v neftegazovyh skvazhinah: pat. RF № 2583867 / M. I. Jepov, V. N. Eremin, A. N. Petrov, V. N. Glinskih ; opubl. 10.05.2016 ; bjul. № 13.
6. Matematicheskoe i fizicheskoe modelirovanie signalov jelektromagnitnogo zonda dlja izuchenija makroanizotropii osadochnyh otlozhenij [Jelektronnyj resurs] / M. I. Jepov, V. N. Glinskih, V. N. Eremin i dr. // Geomodel' 2017: 19-ja konferencija po voprosam geologorazvedki i razrabotki mestorozhdenij nefti i gaza - Gelendzhik, 11-14 sentjabrja 2017 g.: tezisy dokladov. (43809). - URL: http://www.earthdoc.org/publication/publicationdetails/ ?publication=90886 (data obrashhenija - 01.03.2018).
7. Electromagnetic tool for high-resolution logging: theoretical and experimental studies [Electronic source] / M. I. Epov, V. N. Glinskikh, V. N. Eremin et al. // SPE Russian Petroleum Technology Conference. Moscow, October 16-18, 2017. - P. SPE-187904-MS. - URL: https://www.onepetro.org/conference-paper/SPE-187904-MS (data obrashhenija - 01.03.2018).
8. Martakov S. V, Jepov M. I. Prjamye dvumernye zadachi jelektromagnitnogo karotazha // Geologija i geofizika. - 1999. - T. 40, № 2. - S. 249-254.
9. Samarskij A. A. Teorija raznostnyh shem. - M. : Nauka, 1989. - 616 s.
10. Kuznecov Ju. I., Agapitova N. S. Matematicheskie osnovy modelirovanija na JeVM. -Ju.-Sahalinsk : Iz-vo JuSIJePI, 2003. - 214 c.
11. Voevodin V. V., Kuznecov Ju. A. Matricy i vychislenija. - M. : Nauka,1984. - 318 c.
12. Sogabe T., Zhang S. A COCR method for solving complex symmetric linear systems // Journal of Computational and Applied Mathematics. - 2007. - Vol. 199. - P. 297-303.
13. Surodina I. The GPU Solvers for High-Frequency Induction Logging // Numerical Analysis and Its Applications, NAA 2016, Lecture Notes in Computer Science. - Vol. 10187. -P. 640-647. - URL: https://link.springer.com/content/pdf/10.1007 %2F978-3-319-57099-0.pdf
14. Labutun I. B., Surodina I. V. Algorithm for Sparse Approximate Inverse Preconditioners in Conjugate Gradient Method // Reliable Computing (Interval Computations) Journal. - 2013. -URL: http://interval.louisiana.edu/reliable-computing-journal/tables-of-contents.html#Volume_19.
© H. B. CypoduHa, 2018
