научная статья по теме ЧИСЛЕННЫЙ РАСЧЕТ ВОЛНОВЫХ ПРОЦЕССОВ В ТРЕЩИНОВАТЫХ СРЕДАХ НА ГЕКСАЭДРАЛЬНЫХ СЕТКАХ СЕТОЧНО-ХАРАКТЕРИСТИЧЕСКИМ МЕТОДОМ Математика

Текст научной статьи на тему «ЧИСЛЕННЫЙ РАСЧЕТ ВОЛНОВЫХ ПРОЦЕССОВ В ТРЕЩИНОВАТЫХ СРЕДАХ НА ГЕКСАЭДРАЛЬНЫХ СЕТКАХ СЕТОЧНО-ХАРАКТЕРИСТИЧЕСКИМ МЕТОДОМ»

ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 3, с. 512-522

УДК 519.634

ЧИСЛЕННЫЙ РАСЧЕТ ВОЛНОВЫХ ПРОЦЕССОВ В ТРЕЩИНОВАТЫХ СРЕДАХ НА ГЕКСАЭДРАЛЬНЫХ СЕТКАХ СЕТОЧНО-ХАРАКТЕРИСТИЧЕСКИМ МЕТОДОМ1

© 2015 г. В. И. Голубев, И. Б. Петров, Н. И. Хохлов, К. И. Шульц

(141700Долгопрудный, М.о., Институтский пер., 9, МФТИ) e-mail: w.golubev@mail.ru Поступила в редакцию 27.05.2014 г.

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

Ключевые слова: гиперболические уравнения, сеточно-характеристический численный метод, математическое моделирование, трещиноватые среды.

DOI: 10.7868/S0044466915030096

ВВЕДЕНИЕ

По-видимому, прямой метод характеристик был изначально предложен в работе [1]. В дальнейшем, он был обобщен на случай двумерных [2] и трехмерных [3], [4] задач. Основным его недостатком является сгущение расчетных узлов в отдельных областях расчетной области, что приводит к снижению точности расчетов вдали от них. В дальнейшем в работах в [5], [6] был предложен метод обратных характеристик (сеточно-характеристический метод) на фиксированных расчетных сетках. Вышеперечисленные методы успешно применялись для получения гладких решений задач газовой динамики. Для обеспечения возможности корректного расчета разрывных решений были предложены сеточно-характеристические методы с явным выделением разрывов (см., [7], [8]).

Одним из важных свойств численных методов является монотонность (см. [9]). К монотонным схемам относится, например, схема первого порядка аппроксимации (см. [10]). В [11], [12] предложены монотонные гибридные схемы повышенного порядка точности. Общий подход к построению таких схем с использованием пространства неопределенных коэффициентов был предложен в [10], [13]. Также для повышения порядка сходимости без расширения расчетного шаблона в [14] предложена монотонная схема с использованием продолженных систем.

Изначально сеточно-характеристический метод был применен для задач газовой динамики, а к задаче описания динамических процессов в сплошных средах он адаптирован сравнительно недавно. Среди первых работ, посвященных его применению для численного решения многомерных уравнений динамики упругих, а также упругопластических и вязкоупругих, сред можно отметить работы [15], [16]. В работе [17] с использованием сеточно-характеристического метода решались трехмерные волновые задачи теории упругости. Работа [18], по-видимому, является одной из первых отечественных работ, в которой рассматривались волновые задачи геофизики.

1) Исследование выполнено в рамках работы по гранту РНФ № 14-11-00263 и при частичной финансовой поддержке стипендии Президента РФ СП-2548.2013.5.

В [19] сеточно-характеристический метод обобщен на случай численного решения вязкоупругих волновых задач, в [20] — для исследования сложных волновых процессов в средах многослойной структуры. В дальнейшем разработаны методы до пятого порядка точности, как на гексаэдраль-ных (см. [21]), так и на тетраэдральных расчетных сетках для трехмерных задач (см. [22]). В [23] сеточно-характеристический метод на треугольных сетках расширен на случай моделирования геологической среды с явным выделением неоднородностей (в том числе трещин) и измельчения расчетной сетки вблизи них для повышения точности расчета.

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

ОПИСАНИЕ МАТЕМАТИЧЕСКОЙ МОДЕЛИ И ЧИСЛЕННОГО МЕТОДА

Рассмотрим основные уравнения линейной динамической теории упругости, которым подчиняется состояние бесконечно малого объема линейно-упругой среды:

ру\- = V ¡<3ц,

. (1)

= ЦщЬ к1,

где р — плотность среды, — компоненты вектора скорости смещения, а у и бк1 — компоненты тензоров напряжений Коши и деформации, Vу — ковариантная производная по у -й координате. Вид компонент тензора 4-го порядка определяется реологией среды. Для линейно-упругого случая они имеют вид

1 = ^8у5кг + и(8;-к8 А + 8Я5 ук).

В этом соотношении, которое обобщает закон Гука, X и ц — параметры Ламе, а 8у — символ Кронекера.

Первая строка в системе (1) представляет три уравнения движения, вторая — шесть реологических соотношений. Вектор искомых функций, состоящий из 9-ти компонент, имеет вид

и = СТШ СТ12, СТ13, ^22, ^ СТ33Г.

Отметим, что система динамики деформируемого твердого тела (1) может быть записана в матричном виде:

дИ = У Л; , (2)

дг , дх:

у=1 1

где Л у — матрицы размера 9 х 9, (х1, х2, х3) — ортонормированная система координат.

Для ее решения используется сеточно-характеристический метод на криволинейных гексаэд-ральных расчетных сетках из [21]. Исходная система уравнений (2) преобразуется к виду

I=у л у ди, Л у=у ^ Л, , (3)

дг 1=1 у = дх<

где , 2,2, 2,3) — направления расщепления в преобразованном пространстве. Система (3) распадается на три одномерные системы уравнений:

ди = Л,^, у = 1,2,3. (4)

дг 1 д%/

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

ди ,-»-1« ^ ди

= ß -Л ß

д?.

dt 1 1 1 д?

где матрица Q — матрица, составленная из собственных векторов, Л j — диагональная матрица. Для каждого шага расщепления для направления 2, j матрица Л j имеет вид

Л j = diag{ cilj, -e1lpC2lp-c2li,c2li- c2li,0,0,0},

где

c c2 = lj = jwjj = V(w/)2 + (W)2 + (w3)2, wj = V.^..

Для определения коэффициентов w' вычисляется матрица, обратная к матрице Якоби. Отметим, что ее необходимо вычислять в каждой ячейке. Коэффициенты могут быть вычислены аналитически в случае, когда преобразование задано в виде \ = \ (x), либо численно, используя формулы второго порядка точности (wj)m = ((^)m+1 - (^j)m_1)/Axi.

После замены переменных v = Q u каждая из систем (4) распадается на девять независимых скалярных уравнений переноса (индекс j далее опускается, где это возможно) вида

^ + Л-ÖL = 0.

dt d^j

Одномерные уравнения переноса решаются с помощью метода характеристик (см. [8, 24]). После того, как все компоненты v перенесены, восстанавливается само решение:

n+1 n +1

u = fi v .

В программе реализованы схемы 2-4-го порядка точности. В работе расчеты проводились с использованием схемы 4-го порядка точности

v m+1 = v m - tf(Ai - ст(Л 2 - ст(Дз - СТЛ4))),

д=:1(-2vm+2+16 vm+i -16 vm-i+2 ^^

24

Л 2 = 24(-v m+2 +16 v m+1 - 30 v m +16 v m-1 - v m

Дз=:1(2v m+2- 4v m+1+4v m- 2

24

д 4=:1( vm+2- 4v m+1+6vm- 4 vm-1+

24

Кроме того, используется сеточно-характеристический критерий монотонности. В случае положительного значения X критерий монотонности имеет следующий вид

( n n 1 ^ n+1 ^ inn

min jvm vm-1} < vm < max{ vm, vm

m v m-1

Для отрицательных значений А он будет симметричен. В простейшей реализации в случае невыполнения данного критерия монотонности применяется корректировка решения вида

n+1

max

{n n ) n+1 i n n )

Vm^m-} Vm > max {Vm, Vm-1} ,

( n n ) n+1 ^ ♦ ( n n )

mln{Vm,Vm-l), Vm < min{ Vm-1} ,

n+1 ( n n } ^ n+1 ^ ( n n )

Vm , min{Vm-1} < Vm ^ maX{Vm> Vm-1} .

Данный ограничитель сохраняет 4-й порядок в областях, где решение ведет себя достаточно гладко (выполняется характеристический критерий). В случае больших градиентов решения порядок схемы понижается до 3-го.

Приведенный выше алгоритм позволяет находить решение во внутренних точках области. Для этого случая все характеристики, выходящие из расчетного узла на следующем временном слое, остаются лежать в области интегрирования. В граничных узлах ситуация иная. Одна характеристика выходит из расчетной области (см. [25]). В таком случае необходимо использовать граничные условия. Учет граничных условий осуществляется практически одинаково для всех границ. Рассмотрим переход от инварианта Римана к естественным переменным в граничном узле:

П+1 И+1 \ ' П+1 /г\

u = a v = ^Vi r;-, (5)

i

где r-векторы — столбцы матрицы fi"1, v"+l — соответствующие компоненты вектора v"+1. Разобьем сумму (5) на две части. Для первого слагаемого оставим те члены, для которых соответствующие характеристики попадают внутрь расчетной области, обозначим множество индексов i через Xin. У второго оставим те характеристики, которые выходят за границы области интегрирования, т.е. Xout. Тогда равенство (5) можно переписать в виде

n+1 Ч 1 n+1 Ч 1 n+1 Ч 1 n+1 Ч 1 n+1 ъ OUt out in ъ OUt out /f\

u = r = ^ Vn r + ^ Vn r = ^ V< r + fi v = u + fi v , (6)

ieT ieXoU- ieX"'

где fiOUt — прямоугольная матрица, составленная из правых собственных векторов, соответству-

ющих характеристикам, которые выходят за область интегрирования. Вектор V — вектор, составленный из соответствующих компонент инварианта Римана. Вектор и1П может быть найден решением соответствующих уравнений переноса (4). Остается найти вектор V . Граничное условие в линейном случае может быть представле

Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.

Показать целиком