CC BY
74
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: [email protected]
Евгений Игоревич Роменский,
Институт математики им. С.Л. Соболева СОРАН, ул. Коптюга, 4,
630090, г. Новосибирск, Россия E-mail: [email protected]
