О детектировании шероховатых границ

УДК 539.3, 550.3, 517.968.28, 517.956.224

О детектировании шероховатых границ

Е.Б. Сибиряков

Институт нефтегазовой геологии и геофизики им. А.А. Трофимука СО РАН, Новосибирск, 630090, Россия

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

Ключевые слова: метод граничных элементов, упругие стационарные колебания, шероховатые границы, дифракция, медленные волны

On the detection of rough interfaces

E.B. Sibiriakov

Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, 630090, Russia

A new numerical method of solving mixed-type boundary value problems for elastic stationary oscillations has been used, which is based on a boundary element method modification. A full displacement field induced by a two-layered medium with rough and smooth interfaces is calculated in a wide frequency range. It is shown that the signal from an axially symmetric source may contain transverse components after the reflection from the rough interface. These components caused by diffraction only from rough surfaces slowly attenuate with distance from the source and produce slow waves in the real time domain. This qualitative difference in the reflection properties of rough interfaces can be used for their detection.

Keywords: boundary element method, elastic stationary oscillations, rough interfaces, diffraction, slow waves

1. Введение

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

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

Один из первых эффективных инструментов вычисления отражения от негладких поверхностей заключался в использовании лучевого разложения [1]. Однако в случае достаточно быстрого изменения вектора нормали к поверхности это давало очень большие погрешности, чего авторы не скрывали. Аналитические

© Сибиряков Е.Б., 2017

методы [1], использующие метод малых возмущений (метод Рэлея, Кирхгоффа или их модификации), предполагают предельно простую геометрию поверхности раздела. Главное, чтобы вектор нормали к поверхности изменялся не очень быстро. Использование метода конечных элементов для подобных границ раздела затруднительно, поскольку приводит к плохо обусловленным системам, что, в свою очередь, приводит к получению решений с сомнительной достоверностью [2] даже после проведения регуляризационных процедур.

Метод граничных интегральных уравнений предназначен для того, чтобы найти компоненты вектора перемещений и. (либо нагрузок р.), удовлетворяющие как уравнению стационарных колебаний (1), так и граничным условиям:

д и + к и + (Шуи) = 0, к2 = Р^-

ц ц

где X и ц — упругие модули Ламе; р — плотность тела; ю — циклическая частота стационарных колебаний. Решение ищется в виде свертки с некоторым ядром, которое аналитически удовлетворяет уравнению в любой фиксированной точке замкнутого объема х0. Если у нас есть поверхность, на которой нужно выполнить граничные условия, то суммирование проводится в точках этой поверхности х (бегущая точка):

и.

2п

Тензор Мк (х0, х) удовлетворяет системе (1) в фиксированной точке объема х0. Если, например, на поверхности задан вектор перемещений, то фиксированную точку х0 следует устремить на поверхность. В этом случае (2) будет системой линейных интегральных уравнений (т.е. множества х0 и х будут совпадать). Это даст возможность вычислить компоненты вектора потенциала Fk. Если же вектор потенциала уже известен, то это даст возможность в соответствии с (2) вычислить перемещения в любой точке внутри объема, ограниченного поверхностью. Если поверхность является бесконечно удаленной, то на ней потенциал автоматически равен нулю. По этой причине упомянутый выше метод может представлять интерес для подавления ложных отражений.

Классический подход [3] предлагает в трехмерном случае выбрать в качестве М.к (х0, х) тензор фундаментальных решений для пространства:

, ^ / ч т^ / ч 1 (соъ(кг) „

М.к(х0, х) = Гк(х0, х) = -| -+

к2 дх. дхк

ц

со$

где 8.к — символ Кронекера; г — расстояние между фиксированной и бегущей точками; у = ^ц/(Х + 2ц);

х. — проекции вектора, направленного из фиксированной точки в бегущую, на оси декартовой системы координат.

В статике (т.е. при кг << 1) этот же тензор имеет вид

1

Iк (х0, х ) =

2ц(Х + 2ц) дг дг

(Х + ц)—— + (X + 3ц)8й

дх. дхк

г

Его недостатком является то, что тензор нагрузок от него дает, кроме 8-функции, еще и неинтегрируемую особенность. В самом деле, нагрузки от (4), вычисленные в фиксированной точке, имеют вид

Рх0(Гк)(х0, х) =: 1

ц8.к + 3

Х + 2ц дг дг

дх. дхк

д

дп0

1

+ ц

со8

дх. г

1 /- ^

дхк

где п0 — единичная нормаль в фиксированной точке поверхности.

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

2. Формулировка задач

Верхняя свободная поверхность — плоскость ^ =0). В качестве границы раздела использовалась либо плоскость ^ = -0.5), либо шероховатая поверхность ^ = = -0.5 + 0.01(соз(10лх) + зт(10лу))) (рис. 1). Параметр к изменяется от 0.1 до 20, что соответствует эффективной длине 5 волны от 20п до п/10. Упомянутая частотная область охватывает диапазон от квазистатики до лучевой динамики. Нагрузка на верхней поверхности осесимметрична и имеет только нормальную компоненту, равную рп = ехр(-1000г2), г е [0, 1]. Упругие константы верхнего слоя: X = ц = р = 1. Нижний слой имеет такие же параметры, за исключением ц- = 1.2. Условия на границе раздела — непрерывность вектора перемещений и нулевая векторная сумма вектора нагрузок. Использовались цилиндрические координаты с

Рис. 1. Шероховатая поверхность раздела z X

-0.5 + 0.01 X

равномерной сеткой 25x25 на каждой из поверхностей. На верхней поверхности вычислялся вектор перемещений, порожденный нагрузкой и двухслойной средой при различных частотах k.

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

В [5] краевые задачи смешанного типа решались с помощью модифицированного метода граничных интегральных уравнений для статических задач. Модификация заключалась в использовании в качестве ядер соответствующим образом ориентированных статических фундаментальных решений для полупространства. Более сложная конструкция ядер по сравнению с обычными ядрами для безграничного пространства была оправдана лучшей обусловленностью линейных уравнений, соответствующих интегральным уравнениям. Линейные системы для нахождения вектора потенциала, к которым свелись интегральные уравнения, оказались хорошо обусловленными, несмотря на то что такие задачи формально считаются некорректными. Использование модифицированных ядер позволяет существенно расширить класс некорректных задач, которые можно решать без применения регуляризационных процедур. Теперь необходимо обобщить указанный подход и построить гладкие и финитные ядра для случая стационарных колебаний.

3. Построение новых ядер

Можно было бы предположить, что при использовании в качестве ядер фундаментальных решений для

полупространства, можно избежать главного недостатка, связанного с появлением неинтегрируемого ядра, что приводит одновременно и к потере точности, и к повышению числа обусловленности системы. На любой поверхности есть одно выделенное направление — вектор внешней нормали к поверхности, на которой заданы граничные условия. По этой причине в качестве координаты x, = (г, n) выбирается скалярное произведение радиус-вектора, начало которого находится в фиксированной точке объема, а конец — в бегущей точке поверхности, на нормаль в бегущей точке. В зависимости от знака Xj ядрами будут решения для верхнего или нижнего полупространства. Также на поверхности есть два взаимно ортогональных касательных направления eTj, eT. Обозначим

x2 = (r > eт, )> x3 = (r > eт2 )> rT =Vx2 + X2>

x2 = rT cos Ф, x3 = rT sin ф.

Однако тензор Грина для полупространства в цилиндрических координатах [3] в случае стационарных колебаний сохраняет практически все недостатки аналогичного тензора для безграничного пространства. Решение задачи о нормальном ударе по полупространству (задача Лэмба) имеет вид ~ , 2

Un = 1 )_

0 2n^R

(2kr exp(- | x, | v) -

- (2kr2 -k2)exp(- | x, | v ))dkr,

Ur = f J (krrT)-kr-

0 2n^R

(-2v,v exp(- | xj| v,) +

+ (2k2 - k2)exp(- | x, | v ))dkr.

При этом V ^ = -\ кГ - к2, V р кГ -у2 к2; кг — параметр интегрирования; R = 4к2^к2 - к2 -у/к2 - у2к2 --(2кг2 - к2)2 — знаменатель Рэлея; у2 = ц/(А + 2ц); JQ (кггт) — функция Бесселя первого рода. Упомянутые недостатки очевидны. Во-первых, интервал интегрирования бесконечен (потеря точности). Во-вторых, в некоторой точке знаменатель Рэлея обращается в ноль (еще большая потеря точности).

Модернизация ядер производилась в два этапа: замена бесконечного промежутка интегрирования на конечный, затем сдвиг по пространственным частотам. Первый этап в итоге дал существенное (на несколько порядков) повышение точности вычисления ядер. Дело в том, что дельта-функция представляет собой хороший инструмент для аналитического моделирования. Однако при использовании численных методов нагрузку лучше сосредотачивать не в точке, поскольку численное моделирование ограничено пространственными частотами Найквиста (обратно пропорциональными шагу разбиения), а на одной (или нескольких) малых площадках. Это означает, что интегрировать в (7) нужно не до

бесконечности, а до обратного характерного размера разбиения. Другими словами, вместо классической поверхностной дельта-функции при использовании численных методов правильнее применять ее конечный аналог, имеющий ограниченный пространственный спектр (т.е. фактически это дельта-функция, «размазанная» на несколько элементарных ячеек). Заметим, что (7) есть отклик по перемещениям на нормальную нагрузку, представленную в виде (касательные компоненты равны нулю):

8(5) = -—J krJo(krrz )dkr

2п

Сказанное выше о конечном аналоге дельта-функции означает, что в (8) предел интегрирования будет не бесконечность, а некоторое конечное число N. Это означает, что «излишнюю» обусловленность мы «поменяем» на дополнительную точность.

Чтобы избавиться от проблем, связанных с обращением в ноль знаменателя Рэлея, нужно выбрать иное приближение к дельта-функции, поскольку ее представление не единственно. Например,

8( х) = -А1,

дх1 r в то же время

. д cos(kr) 8( х) = --— •

дх1 r

Если первое представление более уместно для статики, то второе — для стационарных колебаний. Каким образом можно исправить (8), чтобы вычислять ядра быстро, точно и достоверно? Ответить на этот вопрос можно с помощью следующей аналогии. Отличие 1/ r от cos(kr)/r в том, что первая формула удовлетворяет уравнению Лапласа, а вторая — Гельмгольца. Параметр уравнения Гельмгольца k во второй формуле примерно одинаковым образом действует на все три измерения. Поскольку в (8) 8-функция поверхностная, т.е. rT заключает в себе две координаты, а не одну, то представляется целесообразным скорректировать (8), учитывая предыдущие замечания о пределе интегрирования, следующим образом:

1 N I-

1 " "

1 1 I

81 (5) = -— J WWk2 + 2k2 )dkr •

2n 0

Отметим, что существенное повышение точности следует из гладкости и финитности ядер. Это дает возможность изменить (7) и получить два из пяти членов матрицы ядер:

N к V

ип =1 Л^згт --

2nyR1(kr, k) х(2v2 exp(- |x11 v1)- (2kr2 + 3k2)x x exp(- | X1I V2))dkr,

N k V

Ur = J J0bRV2 , ) X

0 2n^R1

x (-2v1v3 exp(- | x11 v1) + + (2k2 + 3k2) exp(- |x11 V2))dkr,

V1 =4 k2 + k 2, V2 =4 k2 + k 2 (2-у2),

V3 = *ikl + 2k2,

{k

R1 = 4

- (2к2 + 3к2)2. Важным здесь является то, что знаменатель R1 нигде не обращается в ноль (за исключением случая к = 0 (статика), где нужно раскрыть неопределенность). Кроме того, решение (10) ограничено при х1 = гт = 0. Если число N устремить к бесконечности, то (10) будет стремиться к классическому решению об отклике упругого полупространства на сосредоточенную нормальную нагрузку. При этом оно будет более точно вычисляемым.

Осталось использовать решение для касательного монохромного удара по полупространству, изящно изложенное в [6]. Проведем те же вычисления, только будем искать отклик не на нагрузку (8), а на (9) и используем более удобную систему обозначений. Если х1 > 0 и касательная нагрузка прикладывается в направлении еТ1 (оно же совпадает с направлением оси х2), то вектор перемещений в полупространстве имеет вид:

\кг ехР(- 1 х1| ^ ^

M

N

х 2 ~ J I ~

0 [ 2n^v1

" 0 Jo(V3rT) + J0 (V3rT )er cos ф _ J0(v3rT) еф sin Ф

kr (vj2 +V32)exp(- | X1I V1)) x x

2n^R1

-V3 J0(v3rT)en cosФ

V1J0

-v2 J0 (v3rT)en cosФ

V3 J0

V3[ J0(V3rT) + J0(V3 rT )>ф sin ф_ где er и еФ соответствуют радиальному и азимутальному направлению цилиндрической системы координат, при этом связь с координатами, использовавшимися ранее, определяется соотношениями (6). Теперь найдем все ненулевые компоненты конечного аналога тензора фундаментальных решений для полупространства.

Если перейти к компонентам еп, ег , еф (соотношение между еГт, еф и еТ1, еТ2 определяется (6)), то увидим, что

Мпп =1 J0

к„ V

т" 2

2пцЯ1

х^з ехр(-1х,| V!)-- (VI2 +vз)exp(- | х,| V2))dkт,

I , к v3

Мпт = 1 J0 ^зГ)-—X 0 2л^

х (-2VlV2exp( -1Х11 VI) +

+ (V2 +vз)exp(- |Х11 V2))ёкт,

Мтп = 1J0 ^ )

к„ V

т3

2пцЯ1

х^^2 ехр (-1Х11 V2) -

- (V12 + "^ехр^ 1х^ V 2 ))dkт,

ч к V

Мтт =1J0 (VзTт )-—кг?1—тТ х

х((2vз + vз)exp(-1Х11 Vl)-

- 2Vз2exp(- | Х11 V2))dkr +

+ I

ч к v

Мфф = !-

х((Vl2 + V2 exp (-1Х11 Vl) -

- 2Vз2exp(- | Х11 V2))dkт -

IТ"

Т0 ^з Тт)-;-dkт •

0

Остальные компоненты тензора равны нулю: Мпф = = Мфп = Мфт = Мтф = 0. Таким образом, ядра для решения задач второго рода (на поверхности заданы нагрузки, необходимо найти перемещения) построены. Тензор М к (х0, х) имеет только пять независимых ненулевых компонент в координатах, связанных с направлениями п, гт и ф. Представление ядер в виде конечных аналогов фундаментальных решений приводит к тому, что решение системы (1) естественным образом ищется в виде конечной суммы, а не интегралов (2), т.е. решение ищется в виде

и

Вектор нагрузок можно вычислить аналитически с помощью закона Гука р1 (х0) = ^гкп0к:

р1

Каким образом решения для полупространства (13)-(17) можно использовать для нахождения решений краевых задач на произвольных поверхностях? Если кратко,

то ядра нужно правильно сориентировать. Что это означает? В формулах

ип

Тт 1 Тт

Однако уравнение (20) использовать невозможно. Вместо него следует использовать проекцию в фиксированной точке (а в ней направления нормали и касательных иные, чем в бегущей):

и п0(х0) = Е Мп0п (х0> х)Р'п (х) + + Мп0Т1 (Л0, х)FTl (х) + Мп0т2 (Л0, х)FTз (х). (21) При этом компоненты тензора

следует вычислять как соответствующие проекции вектора, например:

х2 Тт

М,

.0п = Мпп (П0' п) + Мп (п0> О" +

+ Мтп

хз

т.е. все необходимые проекции М к (х0, х) могут быть вычислены как компоненты вектора после применения операции вращения. В то же время необходимые компоненты тензора нагрузок (Рк (х0, х) из (19)) вычисляются как компоненты тензора второго ранга после вращения. Обратим внимание на то, что новые ядра в форме (13)—(17) не зависят ни от вида граничных условий, ни от геометрии поверхности. Из всех величин, определяемых геометрией поверхности, они зависят только от | х1| и тт, точнее от их максимальных значений на поверхности. Поэтому перед решением конкретной упругой задачи (и особенно серии задач) имеет смысл сначала определить максимальные расстояния между точками, которые возможны на этой поверхности (эти расстояния и будут верхней границей для | х11 и тт, нижняя граница — ноль). После этого можно создать библиотеку, которая будет представлять собой эти пять двумерных функций (13)—(17), вычисленных на определенной сетке, а также их частные производные, вплоть до третьих включительно. Таким образом, при

х

т

Рис. 2. Угловая компонента еф вектора перемещений на верхней границе, порожденная шероховатой границей раздела

0.0 0

Рис. 4. Угловая компонента еф вектора перемещений на верхней поверхности, порожденная отражением от шероховатой поверхности

решении конкретной группы задач на фиксированной частоте стационарных колебаний это дает возможность избежать вычисления в каждой точке упомянутых выше интегралов, используя двумерную формулу Тейлора. Например, если известно, что нужно будет решать задачи, в которых максимальное расстояние между точками будет равно 2 (см. пример ниже), то можно положить число N в (13)—(17) равным 200. Вычисление в каждой точке векторов п, еТ1, ет , их скалярных произведений в фиксированной и бегущей точках, а также параметров х1, х2, х3, гТ и элемента поверхности, а также столбца свободных членов (который зависит от граничных условий) происходит несравнимо быстрее, чем вычисление интегралов. Поэтому процесс вычисления матрицы коэффициентов для системы линейных уравнений, предназначенной для нахождения вектора-потенциала, становится приемлемым по времени (см. пример ниже).

4. Результаты

Отражение осесимметричного сигнала от шероховатой поверхности дает угловую компоненту на всех частотах (эта компонента при отражении от гладкой

поверхности есть ноль). Однако на достаточно коротких волнах амплитуда этой компоненты слишком мала по сравнению с поверхностной волной. Угловая компонента представлена на рис. 2 при к = 10. Однако амплитуда этой компоненты уменьшается при удалении от источника. Разность радиальных компонент, порожденных гладкой и шероховатой границами раздела (к = 10), представлена на рис. 3. На высоких частотах шероховатость может быть выделена вычитанием радиальных компонент, полученных при отражениях от гладкой и шероховатой поверхностей соответственно. Эта разность есть квазипериодический сигнал с незначительным амплитудным затуханием.

Угловая компонента еф вектора перемещений на свободной поверхности, порожденного отражением от шероховатой поверхности, представлена на рис. 4 при к = 0.1. Это означает, что если волна достаточно длинная (квазистатика), затухание угловой компоненты становится незначительным и ее амплитуда становится сравнимой с амплитудой поверхностных волн.

Если вернуться во временную область и построить угловой аналог сейсмограммы, т.е. зависимость угловой компоненты вектора перемещений от времени и угла

0.0 0

Рис. 3. Зависимость разности между радиальными компонентами, порожденными отражениями от плоской и шероховатой границ

Рис. 5. Круговой аналог сейсмограммы, т.е. зависимость угловой компоненты вектора перемещений от угла и времени при фиксированном расстоянии до центра источника иф

при фиксированном (единичном) расстоянии от источника иф (1, ф, 1), то можно выделить волны, соответствующие этому перемещению и определить их скорость. Выяснилось, что при отражении от шероховатой границы в некотором диапазоне частот образуются волны (точнее, пара волн) с очень низкой скоростью распространения, причины возникновения которых рассматривались в [7]. Их вполне можно идентифицировать как медленные волны. Их скорость при данных параметрах численных экспериментов примерно в десять раз меньше скорости поперечных волн, также существенно меньше их геометрическое расхождение (рис. 5).

5. Выводы

Шероховатая граница раздела порождает угловую компоненту вектора перемещений на свободной границе при осесимметричном источнике. В случае если граница раздела плоская, эта компонента отсутствует.

Амплитуда этой компоненты растет с ростом эффективной длины волны. Эта компонента затухает существенно слабее, чем остальные компоненты. Это дает возможность выделить ее на большом расстоянии от источника.

Если длина волны достаточно мала (по сравнению с толщиной слоя), то тогда обнаружить шероховатость

можно с помощью вычитания радиальных компонент, порожденных плоской и шероховатой границами раздела.

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

Возможно, что медленные волны порождаются не объемными, а поверхностными дефектами.

Литература

1. Kennett B.L.N. Seismic wave scattering by obstacles on interfaces // Geophys. J. R. Astr. Soc. - 1972. - V. 28. - P. 249-266.

2. Иванов K.M., Винник П.М., Иванов B.H. Численное моделирование

разделительных процессов обработки давлением // Вестник СГФУ. - 2012. - Т. 33. - № 2. - С. 192-199.

3. Новацкий В. Теория упругости. - М.: Мир, 1975. - 872 с.

4. Купрадзе В.Д. Методы потенциала в теории упругости. - М.: Физматгиз, 1963. - 472 с.

5. Сибиряков Б.П., Сибиряков Е.Б. Области локального понижения давлений как вероятные аккумуляторы флюидов в геологических структурах // Геология и геофизика. - 2015. - Т. 56. - № 7. -С.1391-1397.

6. Ziatdinov S.R., Chestnut B.M. Extrinsic components of the Rayleigh wave // Geophys. Quest. - 2005. - V. 38. - P. 46-55.

7. Сибиряков Б.П. Параметрические резонансы и неустойчивость геологических структур // Физ. мезомех. - 2005. - Т. 8. - № 2. -С. 5-10.

Поступила в редакцию 17.08.2017 г.

Сведения об авторе

Сибиряков Егор Борисович, к.ф.-м.н., снс ИНГГ СО РАН, sibiryakoveb@ipgg.sbras.ru