WENO/Рунге-Кутта метод высокой точности для моделирования упругих волн Текст научной статьи по специальности «Физика»

ISSN 2074-1863 Уфимский математический журнал. Том 2. № 1 (2010). С. 59-70.

УДК 519.6

WENO/РУНГЕ-КУТТА метод высокой точности для МОДЕЛИРОВАНИЯ УПРУГИХ ВОЛН

М.Н. ДМИТРИЕВ, Е.И. РОМЕНСКИЙ

Аннотация. В работе представлено применение численного метода WENO/Рунге-Кутта высокого порядка точности для решения уравнений линейной теории упругости, записанной в виде гиперболической системы законов сохранения. Рассматривались методы до 5-го порядка точности по пространству и до 4-го порядка точности по времени. Сравнение результатов расчетов тестовых задач с результатами, полученными широко применяемым в сейсмике методом Вирье второго порядка по пространству и по времени, показывает несомненное преимущество алгоритма WENO/Рунге-Кутта. Рассмотрено также применение метода в случае ограничения расчетной области посредством введения специальным образом построенных поглощающих слоев PML (Perfectly Matched Layer).

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

1. Введение

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

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

M.N. Dmitriev, E.I. Romenski, WENO/RK method for modelling elastic waves.

© Дмитриев М.Н., РомЕнский Е.И. 2010.

Работа поддержана РФФИ (проекты 07-05-00538, 08-05-00265, 09-05-00221), Президиумом РАН (проект № 2).

Поступила 20 октября 2009 г.

пространству, так и по времени. Можно упомянуть такие методы, как дискретный метод Галеркина, ADER метод, WENO методы [4].

В данной работе представлено применение WENO/Рунге-Кутта алгоритма для решения гиперболической системы законов сохранения линейной теории упругости. Рассматривались методы до 5-го порядка точности по пространству и до 4-го порядка точности по времени. Сравнение результатов расчетов тестовых задач с результатами, полученными широко используемым в сейсмике методом Вирье второго порядка по пространству и по времени, показывает несомненное преимущество WENO/Рунге-Кутта методов. Отметим, что разработанные методы могут быть прямо использованы для расчета упругих волн в средах с переменными (в том числе с разрывными) характеристиками среды (плотность, скорости звука).

В работе также рассмотрено применение метода в случае ограничения расчетной области посредством введения специальным образом построенных поглощающих слоев PML (Perfectly Matched Layer), предложенных в работах [7, 8]. Разработанные алгоритмы иллюстрируются серией расчетов тестовых задач, одномерных и двумерных.

2. Уравнения линейной теории упругости

Мы будем рассматривать уравнения линейной теории упругости как гиперболическую систему уравнений первого порядка, которая формулируется в терминах скоростей движения среды Пг и деформаций £ц.

дрщ да, ~dt д£а

dxi

dt

дщ

dxj

0,

+

i,j

дщ

dxi

= 1, 2, 3,

= 0, i,j = 1, 2, 3.

(1)

Первое уравнение здесь выражает закон сохранения импульса, второе уравнение описывает эволюцию тензора малых деформаций при движении среды. Мы будем рассматривать случай изотропной среды, для которой напряжения а^ связаны с деформациями и законом Гука

= А(£11 + £22 + £зз)^у + ■ (2)

Здесь р — плотность среды, А, ^ — постоянные Ламе.

Известно, что система (1) является гиперболической, что позволяет применять для ее решения современные численные методы, разработанные для решения гиперболических систем законов сохранения.

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

Двумерная система уравнений может быть получена из (1) в предположении, что движение вдоль одной оси координат, например х3, отсутствует. Это означает, что и3 = 0, а значит, и а13 = 0, а23 = 0, а33 = 0. При этом система (1) сводится к нижеследующей векторной форме

Ш д¥(и) , ав(и)

~т+ '

+

0.

дх ду

Здесь консервативные переменные и и конвективные потоки Е, О имеют вид ( ри1 \ ( — а11 \ ( — а12 \

(3)

U

PU2

£11

£22

V £12 /

F(U)

“а12

U1

0

G(U)

V u2/2 /

“а22

0

— U2

(4)

\ -u1/2 J

Далее для конструирования численных методов будут использоваться одномерные варианты системы (3)

d\J 9F(U)

дt дх

ш 0G(U) _

dt ду ’

(5)

а для решения задачи Римана (задачи о распаде разрыва) — их эквивалентные формы, записанные в характеристических переменных.

Для одномерной системы, описывающей распространение волн вдоль оси х, такая система имеет вид

д

dt \ 1 А + 2 ¡і

д_ ’ дх

о-i і ) ± ср— ( щ Т

д ( cs \ д ( cs

77Г «2 =F —1а12 ± Cs— U2 Т —(712

Л + 2v

ой

dt

V

dx

V

0.

(6)

Здесь Ср = л/(X + 2/х)/р и с8 = л]р,/р соответственно продольная и поперечная скорости звука. Система для волн вдоль оси г имеет следующий вид:

д_

dt \ z 1 Л + 2/1

д ( cs

— Иі =F —стіг

dt \ v

°22

, д

±С*ву

Л + 2v

022 = 0,

4- 9

± cs —

дУ

щ Т —012 = 0. V

(7)

0

c

p

0

c

p

Системы (6), (7) могут быть применены для решения задачи Римана, которая в случае одномерного движения вдоль оси х имеет нижеследующую формулировку. Пусть при t = 0

11111 _ РРРРР

на оси х заданы начальные данные щ, и^, єЦ, £22, £12 при х < 0 и и^, щ, є{\, є^, при х > 0. Требуется найти решение при t > 0.

Решение задачи Римана является кусочно-постоянным в плоскости (х,^, разделенной характеристиками ^х/^ = —ср, ^х/^ = —с5, ^х/^ = 0, ^х/^ = с5, ^х/^ = ср. Нас инте-

ресует решение и1,и2, є1і, є22, ЄІ2 в области, ограниченной характеристиками ^х/^ = —с5 и ^х/^ = с5, и выражающие его формулы имеют вид:

1 / Е> Г n 1 Cs

И2 — 2 2 /І ^0"12 0"12^ ’

= 2 (<^11 + а1\) + 2 + 2 с-^"^ _ ^

°І2 = \^12 + ^и) + ^(и2~ и2)•

Отметим, что задача Римана для случая, когда среда слева и справа от разрыва имеет разные материальные характеристики (р1, с^, с1 слева и р1, с^,с1 справа), также может

быть решена методом характеристик, и ее решение в области, ограниченной характеристиками ^х/^ = — с1 и ^х/^ = с^, имеет следующий вид:

pRCpUi + pLCpUi + (J и — (Ти

pLcl; + pRcR

о

pRcRuR + pLcLuL + aR - an

іі

PLcf + pR c

Rf->R

s

pLCpPRCp

pLcL + pRcR

a

pR cR

p

+

a

ii

PLcL

(9)

a

pLCgpRcf pLcf + ряс^

a

pR cR

s

+

a

i2

PLcf

+ uR —

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

[их ] = 0, [аи ] = 0, [«2 ] = о, Ы = 0,

где [/] = /к — /ь — скачок функции / при переходе через разрыв.

L

L

3. Ограничение расчетной области

Во многих задачах сейсмики и сейсмоакустики волны могут распространяться в неограниченной области. Для того чтобы ограничить расчетную область, применяются различные методы. Ниже будет описана адаптация к применяемому методу высокого порядка точности подхода, использующего окаймление расчетной области поглощающим слоем. Такой слой должен обеспечивать поглощение без отражения волн, приходящих в него из расчетной области. В данной работе мы остановимся на так называемых идеальносогласованных поглощающих слоях PML (от английского Perfectly Matched Layer). Это специальным образом сконструированный слой, расположенный вдоль границы расчетной области и обеспечивающий затухание решения по мере его распространения. Впервые описание таких слоев было изложено в работе [7] для расчета электромагнитных волновых полей, а применительно к уравнениям упругости — в работе [8].

Сформулируем уравнения, моделирующие распространение волн внутри PML и обеспечивающие затухание решения. Решение системы (3) представляется в виде суммы

U = Ux + Uy,

слагаемые которой являются решениями уравнений

-¡W.-P.* ...

Здесь d(s) > 0 — демпфирующая функция, которая, следуя работе [3], выбиралась следующим образом:

2cp , ( s \ 4

Ф) = -^log(l/i?) (j) ,

Ь 04 ' 7 \Ь,

где 5 € [0, Ь], Ь — ширина поглощающего слоя, Я — константа характеризующая коэффициент отражения. В численных расчетах выбиралось Я = 10-5.

4. Численные методы

В данном параграфе приведено краткое описание конечно-объемного метода ШЕКО в применении к описанным в предыдущем параграфе двумерным уравнениям линейной теории упругости (3).

Предположим, что плоскость ж, у разбита на счетные ячейки

Пц = [хг-1/2, жг+1/2] х [у'-1/2, у'+1/2]

с номерами г, и длинами сторон Аж = ж*+1/2 — 1/2, Ау = у^+1/2 — у^- 1/2. Пространствен-

ная аппроксимация уравнения (3) имеет следующий вид

(Ю^ _ 1+1/2,0 - ^1-1/2,0 ^4,0+1/2 - ^4,0-1/2

¿і Ах Ах

Здесь И*/ — значение решения, отнесенное к ячейке Щ/, а Е*_1/2/ ,Еі+1/2/, 0*,/_і/2,СІ;^+і/2 — значения потоков через грани ячейки. Вычисление потоков на гранях счетных ячеек может быть выполнено различными способами. Мы опишем здесь конечно-объемный ШЕКО алгоритм [6], который использует усредненные значения решения по счетным ячейкам:

1 [хі+1/2 |’уз + 1/2

АхАу

= —¡Г- и (і,х,у)сІхсІу.

"г-1/2 У3 —1/2

Выражения для потоковых членов в (11) получаются при интегрировании уравнения по объему ячейки П/ = [х*-1/2,х*+1/2] X [у/-1/2, У/+1/2]:

1 ГУ3 +1/2

¥і+і/2,о = ^— Е(и(і, ж*+1/2,2/))<І2/,

Ау Уу —1/2 1 №+1/2

0+1/2 = д- С(иЦ,х,Уз+1/2))(1х. (12)

Ах J хг—1/2

Опишем общую процедуру нахождения Еі+1/2/ . Для этого интегралы в формуле (12) аппроксимируются ^-точечной квадратурной формулой

1 М

Ті+1/2,о ~ ^2~Р(и(і,Хі+1/2,Уа))Ка, (13)

а=1

где Ка — веса квадратурной формулы. Далее для нахождения величин иг+1Аа = и(^+1/2, уа) применяется ШЕКО реконструкция из средних значений по ячейкам слева и справа на гранях между ячейками иь*+1/2,а = и(^,ж*+1/2 — 0, уа),

ик¿+1/2,а = и(£, Жг+1/2 + 0, уа). После этого приближенное значение Е*+1/2^- получается из

решения задачи о распаде разрыва в узлах уа:

1 М -

Р‘+1Л> “ ^ Р(ин1/*.«' и?+1/2,„Жа. (14)

у а=1

В данной работе для аппроксимации интеграла использовалась двухточечная квадратура 4-го порядка:

^1+1/2,0 ~ 2^ /2,Уо - + 2^ (^(ж*+1/2’ Уо +

2 V х 2л/3 Ч 2

Для вычисления \]{хг+\/21Уо ^ пРименялась двумерная \VENO реконструкция 5-

го порядка. Для этого сначала применяем одномерную реконструкцию и(ж*+1/2, у&, £) из

средних значений и (ж*, у к, £) слева и справа от грани х = ж*+1/2 для к = ] — 3 качестве шаблона используется взвешенная комбинация трех шаблонов:

и(0) = ^(-иг+2 + 5ит + 2иг),

6

и(1) = ^(-иг_! + 5иг + 2иг+1),

6

и(2) = ^(2и*_2 - 7и*_! + 1Ш*),

6

И1+1/2 = И*+1/2 (х*+1/2 — 0) = ^ шгИ(1). Веса ш имеют следующий вид:

1=0

аг ¿о ^1

= а°= (£ + /?о)2’ “1 = (£ + А)2:

1^аг 1=0

¿2

а2 =

(е + 02)2’

Л - 3 А 3 Л - 1

(¿0 — 1 $1 — —1 (¿2 — •

10 5 10

Индикаторы гладкости шаблонов:

13 1

Ро = ~ 2^*+1 и*+г)2 + — 4и*+1 + иг+2)2,

13

А = ^(и*-1 - 2и* + и*+1)2 + (и*-1 - и*+1)2>

13 1

@2 = ~(и*-2 — 2ТХ(_1 + и,)2 + -(и*_2 — 4и*_1 + и,)2.

Формулы ДЛЯ вычисления и 1 имеют вид

г~2

1

и(0) = -(2иг+2 - 7иг+1 + 1Шг), 6

и(1) = ^(-иг+1 + 5иг + 2иг_!), 6

и(2) = и*_2 + 5и*_! + 21] г),

6

и

я

-1/2

и*-1/2(х*—1/2 + 0) = и

(г)

г=о

Веса ¿0, ¿1 и ¿2 получаются циклической перестановкой:

А 1 Л - 3 А 3 ил — ----, Щ — —, Сь2 — ---------•

0 10 1 5 2 10

•+ 3. В

(16)

(17)

(18)

(19)

После этого проводим одномерную ШЕКО реконструкцию по переменной у, используя найденные выше значения. Формулы для реконструкции в точке у^ имеют вид

и(0> = и, + (30, - №(*, + и)+2)^|,

U^11 = Uj — (—UJ_1 + U<2> = Uj - (3Uj - 4U,., + и;_2)^|,

1=0

Веса d0, di и d2 имеют вид

Формулы для вычисления U + ^/0 имеют вид

(21)

, 210 — л/З , И , 210 + л/З

do —----------------, d\ — —, do —--------------------------• V /

1080 18 1080 U

(23)

и"» = и, - (зи, - 4Ц|+1 + иг+2)^-,

и'1» = ц, - (и^ - и;+1)^|,

и® = и, - (-30, + 4и^_, - и,_2)^-,

и(Й + И)^Ш'иИ

Веса ¿0, ¿1 и ¿2 получаются циклической перестановкой:

210^3 11 210 -уз (24)

1080 ’ 18’ 1080 Шаг интегрирования по времени А£ выбирается из условия Куранта-Фридрихса-Леви:

V 13 13

и Syз — максимальные скорости распространения волн в направлении осей х и у соответственно. СХХ — число Куранта-Фридрихса-Леви, которое выбирается как СХХ < 1 в одномерном случае и СХХ < 1/2 в двумерном случае.

5. ЧИСЛЕННЫЕ ПРИМЕРЫ

В данном параграфе будет проведен сравнительный анализ результатов расчетов, полученных с использованием метода WENO/Рунге-Кутта и других классических методов, в частности широко используемого в задачах сейсмики метода Вирье. Применение численных методов в изучении распространения упругих волн сталкивается с такими трудностями, как необходимость точного вычисления волн высоких частот (10-100 kHz) на

очень больших временах пробега (до нескольких сотен длин волн). Мы продемонстрируем преимущество ШЕКО/Рунге-Кутта методов при решении такого рода задач.

5.1. Одномерные упругие волны. Продемонстрируем теперь преимущество метода ШЕКО/Рунге-Кутта по сравнению с некоторыми другими конечно-разностными методами на серии одномерных задач о распространении упругих волн. Рассмотрим вначале волну, распространяющуюся вправо вдоль оси х. Волна задается в начальных данных следующим образом:

иі(0,х) = е-10(х-2)2, £11 (0,х) = -е-10(х-2)2.

Точное решение на момент времени і, соответствующее этим начальным данным, выражается формулой

и1 (і,х) = е-10(х-Ср 1)2, £11 (і,х) = -е-10(ж-Ср*)2.

На рис. 1 слева приведено сравнение результата расчета методом С.К. Годунова распада разрыва (первый порядок точности)[1] с аналитическим решением на момент времени, соответствующий прохождению волной расстояния около 50 длин волн. Первоначально на длину волны приходилось 40 точек, число Куранта в расчетах было взято 0.7. Полученный численно профиль скорости дает сильное падение амплитуды и размазывание профиля. Можно сделать вывод, что такого рода методы дают неудовлетворительные результаты при решении подобных задач.

На рис. 1 справа приведено сравнение с точным решением численных результатов, полученных схемой Вирье и схемой ШЕКО/Рунге-Кутта пятого порядка по пространству и 3-го и 4-го порядков по времени. На длину волны приходится 20 точек, число Куранта равно 0.7. Видно, что схема Вирье создает существенные осцилляции за волной, а сам фронт волны оказывается несколько запаздывающим по сравнению с точным решением. Методы ШЕКО-5/Рунге-Кутта-3,4 дают гораздо лучшие результаты, но падение амплитуды достигает 20%. Можно тем не менее увидеть, что повышение точности интегрирования по времени уменьшает падение амплитуды волны.

На рис. 2 показаны результаты расчета той же задачи с более мелким пространственным шагом (число точек на длину волны — 40), из которых видно, что метод Вирье и в этом случае дает осциллирующее решение, в то время как ШЕКО-5/Рунге-Кутта-4 дает очень хорошее соответствие точному решению.

Метод ШЕКО/Рунге-Кутта легко может быть обобщен на случай переменных коэффициентов уравнений (различные скорости звука и плотности в среде). Оказывается, что для этого не нужно выделять контактную границу, и сквозной счет обеспечивает нужный порядок точности. На рис. 3 слева приведены результаты расчета задачи о распространении волны в двухслойной среде. Внутри расчетной области х Є [0, 100] при х < 50 параметры среды р =1, с = 1, а при х > 50 параметры среды р =1, с = 2. Форма источника прежняя, в начальных данных было взято 20 точек на длину волны. Видно, что метод Вирье дает неудовлетворительные результаты по сравнению с методом ШЕКО-5/Рунге-Кутта-3.

Только при существенном увеличении числа точек на длину волны метод Вирье дает результат, приближающийся по точности к методу ШЕКО/Рунге-Кутта. На рис. 3 справа приведен тот же расчет, но для числа точек на длину волны 20 для ШЕКО-5/Рунге-Кутта-

3 схемы и 400 точек на длину волны для схемы Вирье.

5.2. Двумерные упругие волны. Приведем теперь пример расчета распространения упругих волн, инициированных точечным источником, и их взаимодействия с поглощаю-тттим слоем. Расчетная область представляет собой прямоугольник (х, у) Є [0, 8] х [0, 8]. В начальный момент времени в центре области задается источник в виде импульса Риккера,

который входит как правая часть в уравнения для £11, е22:

де

ii

du

i

dt

де

22

dt

дх

<9u2

ду

/ (t)^(x - xS), / (t)£(x - xS),

(25)

где f(t) = 1 - 27T2z/2(i - i0)2e-^2(i-io)2, г/ = 30, t0 = f.

Поглощающие слои заданы внутри расчетной области в виде (х, у) Є [0, 8] х [0, 1.5], (х, у) Є [6.5, 8] х [0, 8].

Для расчета во всей области, включая PML, использовался описанный выше метод WENO-5/Рунге-Кутта-4. Уравнения, моделирующие затухание внутри PML, брались в виде (10).

На рис. 4 приведено поле напряжения а11, полученное в результате расчета для различных моментов времени t = 1.2, 1.5, 1.7, 1.9. Видно, что волна уходит за пределы расчетной области, при этом не возникает видимых отражений от границы раздела основной области и PML.

X X

Рис. 1. Искажение импульса. На левом графике - распределение скорости, сравнение численного решения по схеме Годунова с аналитическим решением. Справа - сравнение численного решения по схеме Вирье и ШЕКО/Рунге-Кутта методов с аналитическим решением

Рис. 2. Искажение импульса. Сравнение численного решения по схеме Вирье и ШЕКО/Рунге-Кутта методов с аналитическим решением

■ WENO5/RK4 Virieux

1

50 60

X

Рис. 3. Искажение импульса для двухслойной модели. Распределение скорости, сравнение численного решения по схеме Вирье с ШЕКО-5/Рунге-Кутта-4 методом

6. Заключение

Метод WENO/Рунге-Кутта может эффективно применяться для расчета распространения упругих волн, в том числе в средах с переменными характеристиками упругости. При этом точность результатов, полученных методом WENO/Рунге-Кутта, превосходит классические конечно-разностные методы. Включение в расчет поглощающих PML слоев не представляет трудностей в применении метода. Заметим, что представленный алгоритм легко обобщается на трехмерный случай.

Авторы выражают признательность В.А. Титареву, В.А. Чеверде за плодотворные обсуждения данной работы.

X

0

30

40

70

80

X

t = 1.2

t = 1.5

t = 1.7

t = 1.9

Рис. 4. Снимки волнового поля ац для различных времен

СПИСОК ЛИТЕРАТУРЫ

1. Годунов С.К., Забродин А.В., Иванов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. М.: Наука. 1976.

2. J. Virieux P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method // Geophysics. 1986. Vol. 103. №. 4. P. 889-901.

3. E.H. Saenger, N. Gold, S.A. Shapiro Modeling the propagation of the elastic waves using a modified finite-difference grid // Wave Motion. 2000. Vol. 31. №. 1. P. 77-92.

4. E.F. Toro Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer. 2009.

5. G.S. Jiang, C.W. Shu Efficient implementation of weighted ENO schemes // J. Comput. Phys. 1996. Vol. 126. P. 202-212.

6. V.A. Titarev, E.F. Toro Finite-volume WENO schemes for three-dimensional conservation laws // J. Comput. Phys. 2004. Vol. 201. №. 1. P. 238-260.

7. J.P. Berenger A perfectly matched layer for the absorption of electromagnetic waves // J. Comput. Phys. 1994. Vol. 114. P. 185-200.

8. F. Collino, C. Tsogka Application of the perfectly matched layer absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media // Geophysics. 2001. Vol. 66. P. 294-307.

9. C.W. Shu Total variation diminishing Runge-Kutta schemes // J. Mathematics of Computation. Vol. 67. P. 73-85.

10. D. Drikakis, W. Rider High-resolution methods for incompressible and low-speed flows. SpringerVerlag 2004.

Максим Николаевич Дмитриев,

Новосибирский государственный университет, ул. Пирогова, 2,

630090, г. Новосибирск, Россия E-mail: mnd@ngs.ru

Евгений Игоревич Роменский,

Институт математики им. С.Л. Соболева СОРАН, ул. Коптюга, 4,

630090, г. Новосибирск, Россия E-mail: evrom@math.nsc.ru