CC BY
9
УДК 519.632.4+550.832.7
DOI: 10.33764/2618-981Х-2019-2-3-55-62
БЫСТРЫЕ АЛГОРИТМЫ ТРЕХМЕРНОГО ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ПОКАЗАНИЙ ЗОНДОВ ВИКИЗ И БКЗ, УЧИТЫВАЮЩИЕ НЕРАВНОКОМПОНЕНТОЕ ПОЛЕ НАПРЯЖЕНИЙ В ОКРЕСТНОСТИ СКВАЖИНЫ
Ирина Владимировна Суродина
Институт вычислительной математики и математической геофизики СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Лаврентьева, 6, кандидат физико-математических наук, старший научный сотрудник, тел. (913)903-87-41, e-mail: sur@ommfao1.sscc.ru; Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Коптюга, 3, старший научный сотрудник
Галина Владимировна Нестерова
Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Коптюга, 3, старший научный сотрудник, e-mail: NesterovaGV@ipgg.sbras.ru
В работе рассматриваются быстрые алгоритмы трехмерного математического моделирования показаний зондов бокового каротажного зондирования (БКЗ) и зондов высокочастотного изопараметрического каротажного зондирования (ВИКИЗ), выполненные на основе метода конечных разностей. Приводятся результаты численных расчетов для моделей сред с учетом неравнокомпонентых полей напряжений в окрестности скважины.
Ключевые слова: уравнение Гельмгольца, уравнение Пуассона, каротаж, математическое моделирование, конечные разности.
FAST 3D ALGORITHMS OF VIKIZ AND BKZ LOG NUMERICAL SIMULATION TAKING INTO ACCOUNT THE UNEQUAL COMPONENT NATURAL STRESS FIELD IN THE BOREHOLE ENVIRONMENT
Irina V. Surodina
Institute of the Computational Mathematics and Mathematical Geophysics SB RAS, 6, Prospect Аkademik Lavrentiev St., Novosibirsk, 630090, Russia, Ph. D., Senior Researcher, phone: (913)903-87-41, e-mail: sur@ommfao1.sscc.ru; Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, Russia, 630090, Novosibirsk, Koptyug Prospect 3, Senior Researcher.
Galina V. Nesterova
Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 3, Prospect Аkademik Koptyug St., Novosibirsk, 630090, Russia, Senior Researcher, e-mail: NesterovaGV@ipgg.sbras.ru
The paper deals with fast algorithms for three-dimensional mathematical modeling of the Russian lateral log (BKZ) and data of high-frequency induction isoparametric logging tool (VIKIZ), performed on the basis of the finite difference method. Numerical experiments are presented for formation models taking into account unequal component stress-strein state in the borehole environment.
Key words: Helmholtz equation, Poisson equation, logging, mathematical simulation, finite differences.
Введение
Данная статья является продолжением серии работ авторов [1-3]. За время, прошедшее с момента первой публикации, авторами накоплен значительный опыт в этом направлении. Усовершенствованы алгоритмы и программы, используемые для трехмерного моделирования, что позволило существенно ускорить вычисления.
Ранее рассчитывались модели, где учитывалось напряженно-деформированное состояние пород, окружающих скважину, которые имели две оси симметрии - относительно оси X и оси У и временную зависимость. Пласт рассматривался однородным по глубине, от переменной х распределение удельного электрического сопротивления среды (УЭС) не зависело. Теперь исходная модель усложнилась - появилась зависимость и от глубины (переменной х), что привело к необходимости расчетов для каждого конкретного момента времени профилей сигналов каротажных приборов. Таким образом, количество необходимых расчетов существенно возросло.
Постановка задач и методы решения
При бурении и последующих работах на скважине изменяется напряженно-деформированное состояние пород, ее окружающих [4-6]. В случае неравно-компонентного поля напряжений [7] пористость и проницаемость вблизи скважины будут изменяться в радиальных направлениях по-разному в зависимости от азимутального угла. Как следствие этого, распределения водонасыщенности, минерализации пластового флюида, удельного электрического сопротивления (УЭС) в пласте становятся существенно трехмерными.
Типичный пример модели среды, полученной с учетом неравнокомпо-нентного поля напряжений в окрестности скважины, приведен на рис. 1 и 2.
Рис. 1. Модель распределения удельного электрического сопротивления р( х, у, 2) (Омм), рассчитанная для момента времени 24 часа после бурения
Рис. 2. Распределение удельного электрического сопротивления во втором модельном слое (рис. 1) (сверху) в плоскости (х,у)
Итак, распределение электропроводности с( x, y, z) = 1/ p( x, y, z) считается заданным.
Задача ВИКИЗ
Для амплитуд электрического и магнитного поля справедливы уравнения Максвелла:
rotH = сЕ - iwsE - j
- - (1) rotE = i(jH
Здесь j = j0 = 4л-10-7Гн/м - магнитная проницаемость, ю - циклическая частота, среде, s(x,y, z) - диэлектрическая проницаемость, j - источник стороннего тока (плотность тока), Е - электрическое поле, H - магнитное поле.
В предыдущих работах задача формулировалась для уравнения в частных производных второго порядка относительно вектора аномального электрического поля [1-3]:
rotrotEa + Ea (i(s - c)i(ju = E0(с - с0 - ias)i(ju (2),
где с0( x, y, z) - электропроводность в однородной среде, Ea = (EJ, EJ, EJ)- вектор аномального электрического поля в декартовой системе координат (ось Z совпадает с осью зонда и направлена вниз). В качестве источника, мы рассмат-
риваем магнитный диполь с моментом М = {0,0, М2}, направленным вдоль оси 7.
т Г 7""'0 7~'0 7~'0
Компоненты электрического поля Ех, Еу, Ег в однородной среде описываются следующими выражениями:
E0 =
E0 =-
i». у (i+k0 Re -оr
4tZR 2 RK 07
ia>juMz x \tZR 2 R
(1 + k0 R)e
-k,R
E0 = 0,
здесь Я = ^1 х2 + у2 + г2 , к02 = -га/а0.
Преобразуем уравнение (2) с учетом того факта, что в данном случае выполняется тождество - divE = 0. В результате преобразований получим векторное уравнение Гельмгольца:
-A2Ea + Ea (iws - <j)io>ju = E0( j - j0 - ia>s)io>ju
(3).
В рассматриваемом диапазоне частот (от 0.5 МГц до 14 МГц) в соответствии со скин-эффектом электромагнитное поле экспоненциально затухает с удалением от источника. Это позволяет поставить нулевые граничные условия для вектора Ea вдали от источника.
Используя консервативную конечно-разностную схему для случая разрывных коэффициентов на неравномерной сетке [8] приходим к системе линейных алгебраических уравнений (СЛАУ). Далее, как и в [3], симметризуем полученную СЛАУ. Ранее применялся метод эрмитового разложения для решения СЛАУ с комплексной неэрмитовой матрицей. После успешного применения метода сопряженных ортогональных сопряженных невязок (COCR) [9] для двумерных задач ВИКИЗ было решено и в трехмерной задаче применить этот метод решения, как наиболее быстрый.
В качестве предобуславливателя используем оригинальный, полностью параллельный в реализации на GPU (Graphics Processing Unit) алгоритм аппроксимации обратной матрицы (на основе алгоритма Хотеллинга-Шульца) [10].
Полученная матрица СЛАУ в результате дискретизации уравнения (3) имеет большее диагональное преобладание, чем матрица СЛАУ для уравнения (2); соответственно итерационный метод быстрее сходится.
Таким образом, замена исходного уравнения и метода решения позволила существенно ускорить время расчетов. Программа написана с учетом GPU архитектуры. Каждый зонд рассчитывается отдельно, для расчетов требуется 1 видеокарта. Можно рассчитывать, как на кластере, используя одновременно 5 центральных процессоров и 5 видеокарт, так и на ноутбуке (с приемлемой для математических расчетов видеокарте), последовательным образом считая все
зонды. На кластере можно также легко распараллелить процесс для моделей по времени, используя для каждой временной модели свой процессор, а также можно распараллелить и по каротажному профилю. Это выполняется элементарно и не представляет собой особой технологической сложности. В случае больших расчетов это легко можно реализовать. Тестовые расчеты были проведены на 1 видеокарте с использованием одного центрального процессора.
Например, типичное время расчета для зонда с частотой 14 МГц для одного момента времени (длина диагонали матрицы N=2414475, число точек профиля 240, последовательный перебор) составляет 14 минут на видеокарте К40. Расчеты проводились на НКС-30Т ССКЦ СО РАН. Предыдущий аналогичный вариант расчетов занимал 20 минут.
Задача БКЗ
Постановка прямой задачи моделирования показаний зондов бокового каротажного зондирования и метод ее решения на GPU не претерпели каких-либо существенных изменений по сравнению с [3], поэтому лишь вкратце напомним об этом.
Прямая задача моделирования показаний зондов БКЗ сводится к моделированию поля точечного источника в среде с известным распределением электропроводности к(x, y, z) = 1/р(x, y, z), которое описывается уравнением Пуассона. В цилиндрической системе координат данное уравнение имеет вид
1_d_
r dr
f dUa Л 1 d f dUa Л d f dUa Л
к r- н--- К- + — К
I dr J r2 dr I dф J dz V dz j
1А
r dr
К - К
dU dr
0 Л
+
_d_
r2 dr
(К - к)
dU dф
0
+ •
д_
dz
(К - к)
dU dz
0
(4)
где U = U0 + Ua, U - полный потенциал электрического поля, и0 = 1
4яа0 R
по-
тенциал точечного источника, находящегося в однородной среде на вертикальной оси в точке z=0, иа- аномальный потенциал, I- сила тока, Я = л/т^+г2. При удалении от источника потенциал затухает как 1/ Я, поэтому для функции
= 0. Потребуем также выполнения ус-
r=R
= 0, Ua
=±Z
иа вдали от источников иа ловия периодичности иа\ ^=0 = иа\ ^=2л.
Дискретизация уравнения (3) конечно-разностным методом [8] и последующая его симметризация приводит к системе линейных алгебраических уравнений
Ax = b,
где A - действительная, симметричная, сильно разреженная, положительно определенная матрица. Полученная система решается методом сопряженных градиентов с оригинальным предобуславливателем, полностью параллельным в реализации на GPU. В отличие от алгоритмов для ВИКИЗ, алгоритмы для БКЗ позволяют рассчитывать показания всех зондов одновременно. Что касается расчетов, относящихся к разным временным моделям, они могут быть выполнены как на 1 процессоре (1 видеокарта), так и на кластере с распараллеливанием по каротажному профилю и по моделям. Представленные ниже расчеты были выполнены на 1 процессоре с одной видеокартой. Расчет 240 точек профиля для модели, соответствующей одному моменту времени (длина диагонали матрицы N=1683000), занимает всего 4 минуты на видеокарте К40.
Результаты расчетов
На рис. 3, 4 приведены результаты тестовых расчетов показаний зондов ВИКИЗ и БКЗ, соответствующих модели, приведенной на рис. 1. Цель данных расчетов состояла в том, чтобы проверить правильность работы программ на моделях такого типа и примерно оценить время расчетов. Обсуждение самих результатов расчетов, эффект влияния того или иного параметра, отвечающего за неравнокомпонентность поля напряжений в окрестности скважины - это предмет для отдельной статьи.
12 т
0 ..................................................
-4. -2. 0. 2. 4. 6. 8. 1 . 12. 14. 16. 18. 2 . 22. 24. 26. 28. 3 . 32 . 34. 36.
относительная глубина точки записи, м
Рис. 3. Показания зондов ВИКИЗ для модели рис. 1
О
о с с о о к о ф ф
12
10
л> гч43 п43 г* о^ гч43 п43 к/Ь г43 сх43 гч43 п43 к/Ь г43 о,43 гч43 п43 к/Ь
ъ «V V <о- <Ь' »(V- N$5- ^Ь- гр/ гр- ($>' л$1/ о^»'
относительная глубина точки записи, м
8
6
4
2
0
Рис. 4. Показания зондов БКЗ для модели рис. 1 Заключение
Изучение все более сложных объектов, в частности, пластов с различными геомеханическими свойствами не только в радиальном, но и азимутальном направлениях, требует их трехмерного представления. Соответственно, актуальными становятся и быстрые, эффективные методы расчетов каротажных диаграмм, используемых для изучения объекта. В работе представлен усовершенствованный вариант расчета диаграмм ВИКИЗ, который выполняется в 1.4 раза быстрее предыдущего варианта. Подготовлены программы БКЗ и ВИКИЗ для массовых расчетов.
Работа выполнена при поддержке проекта НИР № 0315-2016-0005 «Методы создания, исследования и идентификации математических моделей с помощью суперкомпьютеров», проекта ФНИ № 0331-2019-0015 «Реалистичные теоретические модели и программно-методическое обеспечение геоэлектрики гетерогенных геологических сред», проекта ИСГЗ ФАНО 0331-2016-0034 "Скважин-ная геофизика в электропроводящих анизотропных диспергирующих средах на основе высокопроизводительных решений трехмерных задач, высокоточных данных каротажа и лабораторных исследований керна", программы 1Х.128.3. «Реалистичные теоретические модели и программно-методическое обеспечение магнито-, электродинамики гетерогенных геологических сред».
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Суродина И.В., Нестерова Г.В. Трехмерное численное моделирование показаний зондов ВИКИЗ и БКЗ на графических процессорах // 16-я научно-практическая конференция по вопросам геологоразведки и разработки месторождений нефти и газа «Геомодель 2014», 8-11 сентября 2014 г., Геленджик. - 4 с. DOI: 10.3997/2214-4609.20142232, http://earthdoc.org/publication/publicationdetails/?publication=77926.
2. Моделирование влияния неравнокомпонентного поля напряжений в окрестности скважины на диаграммы ВИКИЗ и БКЗ [Электронная публикация] / Нестерова Г.В., Ельцов И.Н., Назаров Л.А., Назарова Л.А, Суродина И.В. // Тезисы конференции «Геомодель-2014: 16-я научно-практическая конференция по вопросам геологоразведки и разработки месторождений нефти и газа», (г. Геленджик, Россия, 8-11 сентября 2014 г.) - 4 с. DOI: 10.3997/22144609.20142233, http://earthdoc.org/publication/publicationdetails/?publication=77927.
3. Суродина, И.В, Нестерова, Г.В. Моделирование показаний зондов ВИКИЗ и БКЗ на графических процессорах / Петрофизика сложных коллекторов: проблемы и перспективы 2015. - Сборник статей EAGE. - C. 85-94.
4. Интерпретация геофизических измерений в скважинах с учетом гидродинамических и геомеханических процессов в зоне проникновения / Ельцов И.Н., Назаров Л. А., Назарова Л. А., Нестерова Г.В., Эпов М.И. // ДАН. - 2012. - Т. 445. - № 6. -С. 671-674.
5. Скважинная геоэлектрика нефтегазовых пластов, разбуриваемых на репрессии давления в неравнокомпонентном поле напряжений / Ельцов И.Н., Назарова Л.А., Назаров Л. А., Нестерова Г.В., Соболев А.Ю., Эпов М.И.// Геология и геофизика, 2014, № 5-6. - C. 979-990
6. Эволюция геомеханических и электрогидродинамических полей в массиве горных пород при бурении глубоких скважин / Назарова Л.А., Назарова Л.А., Эпов М.И., Ельцов И.Н. // ФТПРПИ. -2013. -№ 5. - С. 37-49.
7. Hickman S.N., Healy J.H., Zoback M.D. In Situ Stress, Natural Fracture Distribution, and Borehole Elongation in the Auburn Geothermal Well, Auburn, New York //Journal of geophysical research. - 1985. - Vol. 90. - № B7. - P. 5497-5512.
8. Самарский А.А. Теория разностных схем. М.: Наука, 1989. - 655 с.
9. 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. 640647. - URL: https://link.springer.com/content/pdf/10.1007%2F978-3-319-57099-0.pdf.
10. Labutun I.B., Surodina I.V. Algorithm for Sparse Approximate Inverse Preconditioners in Conjugate Gradient Method // Reliable Computing (Interval Computations) Journal http://interval.louisiana.edu/reliable-computing-journal/tables-of-contents.html#Volume_18.
© И. В. Суродина, Г. В. Нестерова, 2019
