CC BY
75
УДК 550.834
Е.В. Лысь, В.В. Лисица, Г.В. Решетова, В.А. Чеверда ИНГГ СО РАН, Новосибирск
КОНЕЧНО-РАЗНОСТНОЕ МОДЕЛИРОВАНИЕ АКУСТИЧЕСКОГО КАРОТАЖА В ТРЕХМЕРНЫХ НЕОДНОРОДНЫХ ТРАНСВЕРСАЛЬНО-ИЗОТРОПНЫХ СРЕДАХ С ПОГЛОЩЕНИЕМ
E.V. Lys, V.V. Lisitsa, G.V. Reshetova, V.A. Tcheverda Trofimuk Institute of Petroleum Geology and Geophysics SB RAS Koptyug, 3 , Novosibirsk, 630090, Russian Federation
FINITE-DIFFERENCE SIMULATION OF ACOUSTIC LOG IN 3D HETEROGENEOUS VTI WITH ATTENUATION
On the base of the application of a finite-difference approximation of an initial boundary value problem for elastic wave equations (velocity/stress formulation), numerical method and its algorithmic implementation have been developed in order to perform a computer simulation of sonic logging. The very general statement is dealt with - surrounding medium allowed to be 3D heterogeneous Vertical Transversely Isotropic (VTI) with attenuation and arbitrary source position at any point inside or outside the well. To provide the most precise description of the sharpest interface of the problem - the interface of the well, the problem is studied in cylindrical coordinates with z axis directed along the well.
In order to truncate area of computations two approaches are used: classical version of Perfectly Matched Layer (PML) and extension on the base of optimal grid. Implementation of parallel computations is done via Domain Decomposition. Data exchange between Processor Units is performed with the help of Message Passing Interface (MPI) library. Results of numerical experiments for VTI media are presented and discussed.
Введение
Ключевое предназначение метода акустического каротажа это детальное определение структуры и механических свойств пород в околоскважинной зоне, посредством измерения и изучения волнового поля создаваемого источником сейсмоакустических волн, расположенным в скважине. Эта задача в упрощённой постановке (скважина в однородной изотропной среде) впервые была исследована в работе [2]. Впоследствии множество авторов внесло свой вклад в изучение этой проблемы [5, 6]. Однако вплоть до настоящего времени отсутствует детальное понимание особенностей процессов распространения акустических волновых полей для реалистичных неоднородных трехмерных сред с анизотропией и поглощением. Поэтому численное моделирование является, по нашему мнению, единственной возможностью изучения волнового поля, возникающего при выполнении акустического каротажа. В работе представлена модификация ранее разработанного метода конечно-разностного моделирования акустического каротажа для изотропной вязкоупругой среды [6], позволяющая проводить расчёты для трансверсально-изотропной упругой среды с поглощением. Благодаря специальному виду тензора упругих модулей в случае VTI среды
можно, как и в изотропном случае применить для построения схемы сдвинутую сетку Верьё [8] и как следствие нет необходимости менять основу разработанного ранее алгоритма, включая измельчение сетки по радиусу и азимуту. Главные изменения коснулись двух областей: описания поглощения и ограничения расчётной области.
Поглощение в анизотропной среде
Для введения поглощения предлагается естественное обобщение подхода, представленного в [9].
Для начала напомним, что в изотропной среде поглощение вводится независимо для P и S волн и определяется добротностями ЯР (т) и ^ (т). В
анизотропной среде ввести поглощение независимо для каждой из трёх волн qP, qSV и qSH не представляется возможным. Рассмотрим тензор, связывающий напряжения с деформациями вязкоупругой среды в частотной области:
Сук1(т) = СцЫ (0 Л1 + Яе $цк1 (т) + 11т $уы (т Л
Наличие вещественной и мнимой частей у компонент тензора влечет за собой наличие вещественной и мнимой частей у фазовых скоростей, то есть:
V2 ( т,п ) = V2 ( 0, п ;[1 + Яе а ](т,п) +1т а ](т п )\ ] = \,2,3.
Будем считать, что мы можем измерить эти комплексные скорости для некоторых направлений на каждой временной частоте. Тогда задача отыскания компонент тензора ставится следующим образом:
Восстановить параметры ¡Зт (т) по комплексным скоростям
(распространения и затухания) измеренным по ряду направлений:
а1 (т, пт), т = 1,..., М
Как можно видеть, этот подход обеспечивает одинаковый тип симметрии (тип анизотропии) и для скоростей и для поглощения.
Мы разработали и применили этот подход для УЛ сред со слабой анизотропией [7] и доказали, что он приводит к классической постановке задачи теории возмущений применительно к симметричной проблеме собственных значений и требуется как минимум пять независимых измерений для каждой временной частоты: четыре для qP или qSV и одно для qSH. Затем необходимо построить рациональную аппроксимацию функции добротности посредством применения обобщенной стандартной линейной модели твёрдого тела, как это сделано в [3].
Ограничение расчетной области
Поскольку метод акустического каротажа в первую очередь нацелен на определение структуры околоскважинного пространства и не ставит задач определения сейсмических объектов на радиальном удалении большем нескольких длин волн, естественным будет ограничить область вычислений одним - двумя метрами от скважины. В вертикальном направлении размер расчетной области определяется размерами акустического прибора и не превышает, как правило, десяти - пятнадцати метров. Таким образом, целевая область имеет форму цилиндра, вытянутого по вертикали. Условия на
границах целевой области должны обеспечивать минимальный уровень артефактов с тем, чтобы не исказить информативные волновые поля. Наиболее распространенным способом постановки таких условий на границах расчетной области является метод идеально согласованного слоя (PML, от английского Perfectly Matched Layer) [6], отлично зарекомендовавший себя в изотропных средах. Но для некоторых анизотропных сред он ведет к возникновению неустойчивых решений, допускающих экспоненциальный рост. Критерий устойчивости PML приведён в работе [1]:
PML в направлении r устойчив, если проекции векторов групповой и фазовой скоростей на это направление имеют один знак (рис. 1). Если этот критерий удовлетворяется, то можно применять стандартный PML. В противном случае, область вычислений расширяется на оптимальных сетках [4]. Основным преимуществом таких сеток является возможность использования весьма крупных шагов (3 точки на длину волны) и обеспечение низкого уровня паразитных отражений на границе между регулярной и оптимальной сетками.
Рис. 1. Критерий устойчивости PML
Численные эксперименты
Первый эксперимент демонстрирует неустойчивость PML для среды, состоящей из скважины заполненной жидкостью и трех горизонтальных слоев с толщинами 1 м, 2 м и 1 м (сверху вниз):
1. isotropic р = 2200 kg/ m3, VP = 3400m / s, Vs = 2500 m / s , s = 0 , y = 0 и S* = 0;
2. VTI p = 2500 kg/ m3, V^ = 4400 m / s , VqSV = 2500 m / s, s = 0.091, y = 0.046 и S'* = 0.688;
3. VTI p = 1800kg/m3, VqP = 3500 m/s , VqSV = 2400m /s, e = 0.215 , y = 0.28 и S'* = 0.359.
Индикатрисы фазовых и групповых скоростей в среде 2 представлены на рис. 2. На рис. 3 представлены произведения радиальных и вертикальных компонент фазовой и групповой скоростей. Как видно, для волн qSV существуют углы, при которых произведения и радиальных и вертикальных компонент отрицательны. Следовательно, для этой среды PML по обоим
направлениям неустойчив. На рис. 4а представлены результаты численного моделирования с использованием РМЬ для ограничения расчетной области, а на рис. 4Ь для этих целей использовано расширение расчетной области с помощью оптимальных сеток.
a) b)
Рис. 2. Медленность (а) и групповая скорость (b) для VTI среды: qP (син.) и
qSV (крас.)
a) b)
Рис. 3. Признак устойчивости PML для VTI среды: qP (син.) и qSV (красн.)
Т=0.53 Т=1.06 Т=1.53 Т=2.00 Т=0.53 Т=106 Т=1.53 Т=2.00
а) Ь)
Рис. 4. Численное моделирование (а) РМЬ (100 точек), (Ь) оптимальное
расширение (150 точек)
На рис. 5 представлен график зависимости относительной ошибки на оси г = 0 от числа точек используемых при расширении расчетной области. Как видно, относительная ошибка при использовании 150 точек в расширении составляет 0.1 %.
100 150 200 250 300 350 400
Рис. 5. Относительная ошибка Рис. 6. Фазовые скорости в модели
от числа точек в расширении для разделения сдвиговых волн
Основной целью следующего эксперимента является наблюдение расщепления поперечных головных волн на стенке скважины. Для этого была выбрана модель среды, содержащая скважину диаметром в один метр, заполненную жидкостью и помещённую в упругую УТ1 среду со следующими параметрами:
р = 2500kg/ш\УчР = 3928т/я, = 2055т/ я, е = 0.334, у = 0.575, 8* = 0.73
Для неаксиально расположенного источника падающая Р волна, генерирует головную qP волну и две головные д8¥ (медленная) и д8И (быстрая). Источник располагался на удалении 0.45 метра от оси скважины. На рис. 7 представлены сейсмограммы записанные приемниками, расположенными по окружности на стенке скважины на вертикальном удалении от источника равном 1 метру. Первая трасса соответствует приемнику расположенному наиболее близко к источнику. Сравнивая сейсмограммы Рис. 7а (УТ1) и рис. 7Ь (изотропия), можно достоверно выделить быструю qSИ волну.
Рис. 7. Радиальные смещения на стенке скважины. (a) VTI. Числами
обозначены первые вступления:
1 - головная P-волна: 2 - головная qSV, 3 - головная qSH (?). (b) Изотропия. Звёздочками обозначены первые вступления P- и S-головных волн
Благодарности
Работа была выполнена совместно с Научно-исследовательским центром компании Schlumberger в г. Москве и частично при поддержке грантов РФФИ 06-05-64748, 07-05-00538 и 08-05-00265.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Becache, E. Stability of Perfectly Matched Layers, Group Velocities and Anisotropic Waves/ Becache, E. Fauqueux, S. Joly P. // INRIA, Rapport de recherche n° 4304, Novembre 2001, 35 p.
2. Biot, M.A. Propagation of elastic waves in a cylindrical bore containing a fluid/ Biot, M.A. // Journal of Applied Physics, 1952, 23, 997-1005.
3. Blanch, J.O. Modeling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique / Blanch J.O., Robertsson J.O.A., Symes W.W. // Geophysics, 1995, 60(1), 176-184.
4. Lisitsa, V. 2005. Optimal grids for numerical solution of a wave equation in heterogeneous media/ Lisitsa, V.// Siberian journal of numerical mathematics, 2005, 8(3), pp. 219-229.
5. Qing, H. L. A 3D cylindrical PML/FDTD method for elastic waves in fluid-filled pressurized boreholes in triaxially stressed formations/ Qing Huo Liu, Sinha Bikash K.// Geohysics, 2003, 68(5), 350-357.
6. Kostin, V. 3D Synthetic Acoustic Log for Viscoelastic Media: Finite-Difference Approach/ Kostin V., Pissarenko D., Reshetova G., Tcheverda V. // Extended abstracts of 69th EAGE Conference and Technical Exposition, London, 11-14 June 2007, P096.
7. Thomsen, L. Weak elastic anisotropy/ L.Thomsen // Geophysics, 1986, 51(10), 19541966.
8. Virieux, J. Velocity-stress finite-difference method / Virieux, J.// Geophysics, 1986, 51, 889-901
9. Zhu, Y. Physical modeling and analysis of P-wave attenuation anisotropy in transversely isotropic media/ Y.Zhu, I.Tsvankin, P.Dewangan, K. van Wijk // Geophysics, 72(1), D1 -D7, 2007.
© Е.В. Лысь, В.В. Лисица, Г.В. Решетова, В.А. Чеверда, 2008
