ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 10, с. 1762-1772
УДК 519.634
ПРИМЕНЕНИЕ СЕТОЧНО-ХАРАКТЕРИСТИЧЕСКОГО МЕТОДА НА НЕСТРУКТУРИРОВАННЫХ ТЕТРАЭДРАЛЬНЫХ СЕТКАХ В РЕШЕНИИ ПРЯМЫХ ЗАДАЧ СЕЙСМОРАЗВЕДКИ ТРЕЩИНОВАТЫХ ПЛАСТОВ1)
© 2015 г. В. А. Бирюков, М. В. Муратов, И. Б. Петров, А. В. Санников, А. В. Фаворская
(141700Долгопрудный, М.о., Институтский пер., 9, МФТИ) e-mail:petrov@mipt.ru; donxenapo@gmail.com Поступила в редакцию 22.10.2014 г.
Переработанный вариант 28.04.2015 г.
Исследуется математическое моделирование сейсмических откликов от трещиноватых геологических пластов с использованием сеточно-характеристического метода на неструктурированных тетраэдральных сетках с применением высокопроизводительных вычислительных систем. Используемый метод был разработан для расчета сложных пространственных динамических процессов в сложных гетерогенных средах и отличается точной постановкой контактных условий. Это делает его пригодным для моделирования задач сейсморазведки в том числе сейсморазведки областей с большим числом неоднородностей, какими являются трещиноватые структуры. Применение неструктурированных тетраэдральных сеток делает возможным задание геологических трещин различной формы и пространственной ориентации, что позволяет решать задачи в постановке, наиболее приближенной к реальной. Использование вычислительного кластера позволяет повысить точность расчета, оптимизировав его продолжительность. Библ. 19. Фиг. 10.
Ключевые слова: сеточно-характеристический метод, неструктурированные тетраэдральные сетки, сейсморазведка, трещиноватые среды, сейсмические волны, численное моделирование, высокопроизводительные вычислительные системы, параллельные алгоритмы.
DOI: 10.7868/S0044466915100075
ВВЕДЕНИЕ
Сейчас сейсморазведка является одним из наиболее распространенных способов исследования пород перед непосредственным бурением скважин (см. [1]). Проводимые исследования направлены на определение структуры геологических пластов, а так же возможного расположения месторождений углеводородов. Проведение численных экспериментов позволяет оптимизировать процесс интерпретации данных сейсморазведки, что способствует уменьшению затрат при добыче нефти. При выполнении численного моделирования таких задач необходимо уметь рассчитывать распространение сейсмических волн в геологических средах с большим числом неоднородностей, не имеющих каких-либо строгих закономерностей в расположении и различающихся по форме, размерам, ориентации и физическим свойствам.
Для того чтобы наиболее точно описать такую сложную гетерогенную среду, оптимальным является использование неструктурированных тетраэдральных сеток. Они позволяют задавать неоднородности (трещины) любой формы и пространственной ориентации.
Система уравнений математической модели состояния сплошной линейно-упругой среды (см. [2], [3]) является гиперболической, поэтому оптимальным является применение именно се-точно-характеристического метода (см. [4]—[11]) с использованием интерполяции высоких порядков (см. [12]), что позволяет добиться наибольшей точности расчета волновых процессов.
Примеры использования сеточно-характеристического метода с интерполяцией высоких порядков на тетраэдральных сетках можно найти в [11].
1) Работа выполнена при финансовой поддержке гранта РНФ № 14-11-00263.
1762
Расчет трехмерных задач сейсморазведки требует обработки большого объема данных. Поэтому необходимо применение высокопроизводительных вычислительных систем. В работе использовались ресурсы вычислительного кластера с распределенной памятью, на который были распараллелены разработанные алгоритмы для оптимального использования ресурсов.
Особенностью используемого в работе подхода является непосредственное задание неодно-родностей в области интегрирования, в чем его выгодное отличие от распространенных моделей эффективных сред (см. [13]—[16]). Он дает более детальные картины волнового отклика и позволяет наблюдать качественно новые эффекты. Трещины задаются в виде контактных или граничных условий, наиболее точно приближающих к модели с заданием физических свойств насыщающего трещины флюида.
В работе рассматривается моделирование откликов от систем однонаправленных трещин, высота которых сопоставима с длиной падающего волнового фронта. Детальное изучение таких геологических структур является одной из наиболее важных задач современной сейсморазведки.
1. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ
Согласно [2], состояние линейно-упругой среды описывается системой уравнений
рдV = (V- О )т, (1)
д(а =Х (V-V) + ц (V® V + (V® у)т). (2)
Уравнение (1) является локальным уравнением движения, где р — плотность материала, V — скорость движения, а — тензор напряжений Коши, являющийся симметричным в силу закона парности касательных напряжений. Уравнение (2) получается из закона Гука путем дифференцирования по времени, X, ц — параметры Ляме, определяющие свойства материала.
В работе используются следующие обозначения: д ?а = да/д? — частная производная поля а
по а ® Ъ — тензорное произведение векторов а и Ъ , (а ® Ъ) = а Ъ1; I — единичный тензор второго ранга.
2. ЧИСЛЕННЫЙ МЕТОД
Для численного решения системы (1), (2) используется сеточно-характеристический метод на тетраэдральных сетках (см. [11]), позволяющий строить корректные численные алгоритмы для расчета граничных точек и точек, лежащих на поверхностях раздела двух сред с разными параметрами Ляме и (или) плотностями.
На каждом шаге интегрирования по времени выбираются три произвольных направления, образующих базис, что обеспечивает изотропность метода, и вводится новая система координат
2,2, 23). В ней систему (1), (2) можно представить в виде
д ? q + А хд ^ + А 2д ^ q + А зд ^ q = 0, (3)
где под вектором q понимается вектор, составленный из трех компонент скорости и шести компонент симметричного тензора напряжений
= {{ь ^ ^ ^ ^ ^12}т .
Для каждой из трех систем вида
5 ? q + А15 ^ = 0 (4)
справедливо точное выражение
I
q((1, (2, (з, ? + т) = X( - ^2, (з, ?), (5)
(=1
где Х( — некие матрицы, выражающиеся через компоненты матрицы А1, о1 — собственные значения матрицы А1, т — шаг интегрирования по времени.
Собственные значения всех трех матриц выражаются через плотность и коэффициенты Ляме следующим образом:
X + 2ц
1/2
Х + 2^ И2
,1/2 / \1/2 И] (И
1/2
,0,0,0}
р у I р ; ^р; фу Фу
В (5) действие матриц X, на вектор неизвестных q можно записать в виде
Х1,2
(п • V )п + -1 (00 - О |)п
С1Р
+р (п • V)(( - Сз)N00 + Сз1) + (N00 - О) (2^00 + XI)
А + 2|а
X
3,4
V V
- ^^5,6
С _С_
(6)
V - (п • V) п + -1 (С • п - (N00 - СТ)п) С2Р
+с2р(п ® V + V ® п - 2(п • V)) + п ® (С • п) + (С • п) ® п - 2(С - N00 Для матриц XI выполняется свойство
£ X, = I - £ X,.
(7)
(8)
В (6)—(8) через п обозначен единичный вектор вдоль направления для матрицы А1, а N 00 есть тензор
N00 = п ® п,
А ^ В есть скаляр:
3 3
А -В = ££ Ав.
I=1 ]=1
Используя в (5) интерполяцию высокого порядка и последовательно применяя для каждого из направлений £,2, Ъз3 формулы, аналогичные (5) и соответствующие системе, аналогичной (4), получаем способ нахождения решения на следующем временном слое. В программном комплексе предусмотрено использование интерполяции (см. [12]) от первого порядка до пятого включительно, что позволяет получать высокую точность численного интегрирования решения по пространству. Также применение матриц XI реализовано с помощью двух операторов, что позволяет сократить количество интерполяций для каждой точки и каждого направления с девяти до шести.
С
3. ГРАНИЧНЫЕ И КОНТАКТНЫЕ КОРРЕКТОРЫ
Используемый метод позволяет применять наиболее корректные вычислительные алгоритмы на границах и контактных границах области интегрирования.
Пусть в матричном виде граничное условие записывается в виде
^^2, %ъ, * + т) = й, (7)
где q £,2,2,3, * + т) — значения компонент скорости и тензора напряжений на следующем шаге интегрирования в граничной точке.
Согласно (6) для каждой матрицы А у имеется три нулевых собственных значения, три положительных и три отрицательных. Пусть для определенности вдоль направления характеристики, соответствующие отрицательным собственным значениям матрицы А1 , выходят за пределы области интегрирования.
Тогда на этапе расчета внутренних точек в соответствии с (5) будет вычислено
qln (,%2,* + т) = £X$( - ст,%2,*).
с >0
Матрица Q*'out составляется из собственных векторов, соответствующих отрицательным собственным значениям.
Действие корректора в граничной точке совершается по формуле
q ((, ^з, t + т) = Fqш t + т) + Фа, (8)
причем условие (7) выполняется с таким же порядком, как и порядок интерполяции.
В формуле (8) матрица (D Q*'out) 1 находится так, чтобы выполнялось условие
(D fi*'out )-1D fi*'out = I, а матрицы Ф и F вычисляются по формулам
Ф = Q *'out(DQ*'out)-1, F = I-ФD.
Для решения различных задач используются граничные условия заданной внешней силы, заданной скорости границы, смешанные граничные условия и неотражающие граничные условия, основанные на равенстве нулю выходящих характеристик. Для неотражающих граничных условий выражение (7) примет вид
vout /r
ЦТ / r q (t + т, гГр ) = 0.
4:
4. УСЛОВИЯ НА ТРЕЩИНЕ
Оказалось наиболее оптимальным с точки зрения производительности расчета реальных не-однородностей (см. [16]) задавать трещины в виде контактных или граничных условий. Были разработаны несколько моделей трещин, задаваемых различными условиями. Трещина считается бесконечно тонкой, между створками трещины находится флюид (нефть, сжиженный газ или вода).
Сначала рассмотрим два предельных случая трещины: 1) трещина с разведенными створками, между которыми находится флюид; 2) полностью слипшаяся трещина, флюид отсутствует.
Для первого случая использовалось контактное условие скольжения. Это смешанное контактное условие
(1п 1п / 1п 1п\ \
РСаУа +Р4СИУь - (Оа - Оь ) • р) • Р
Р РаСа1 + РЬСЬ1
где Ур и fт = 0 подставляется в контактные условия для одной створки трещины с внешней нормалью р, а -Ур и fт = 0 — для другой.
Во втором случае применялось условие полного слипания:
V = -1-(
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.