ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 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 рублей.