Численное моделирование упругих волн разрывным методом Галеркина в неоднородных средах Текст научной статьи по специальности «Математика»

УДК 519.63

М. В. Васильева, У. С. Гаврильева

Численное моделирование упругих волн разрывным методом Галеркина в неоднородных средах

СВФУ им. М.К. Аммосова, г. Якутск, Россия

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

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

Работа выполнена при поддержке проекта РФФИ №17-01-00732a и мега-гранта Правительства РФ №14.Y26.31.0013.

M. V. Vasilyeva, U. S. Gavrilieva

Numerical Simulation of Elastic Waves by the Discontinuous Galerkin Method in Heterogeneous Media

M.K. Ammosov North-Eastern Federal University, Yakutsk, Russia

Abstract. The paper deals with the distribution of elastic waves in an heterogeneous medium. The mathematical model is described by a hyperbolic second-order equation for displacements. For the time approximation, an explicit difference scheme is classically used. Approximation of the equations will be carried out using the Galerkin discontinuous method. This method of approximation makes it possible to obtain a block-diagonal matrix of masses and, consequently, it is easy to convert when constructing an effective computational realization. The results of the numerical solution for the two-dimensional problem

ВАСИЛьЕВА Мария Васильевна - доцент-исследователь научно-исследовательской кафедры вычислительных технологий ИМИ СВФУ им. М.К. Аммосова.

E-mail: [email protected]

VASILYEVA Maria Vasilyevna - Associate Professor-Researcher of the Research Department of Computational Technologies, Institute of Mathematics and Informatics, M.K. Ammosov North-Eastern Federal University.

ГАВРИЛьЕВА Уйгулаана Семеновна - магистрант научно-исследовательской кафедры вычислительных технологий ИМИ СВФУ им. М.К. Аммосова.

E-mail: [email protected]

GAVRILIEVA Uygulaana Semenovna - Master's Student of the Research Department of Computational Technologies, Institute of Mathematics and Informatics, M.K. Ammosov North-Eastern Federal University.

are presented for three model problems with heterogeneous soil properties and also taking into account the presence of cracks.

Keywords: elasticity, wave equation, discontinuous Galerkin method, finite element method, mathematical modeling, cracks, fractured medium.

This work was supported by the RFBR project No. 17-01-00732a and the mega-grant of the Government of the Russian Federation No. 14.Y26.31.0013.

Введение

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

При рассмотрении многомерных волновых уравнений в неоднородных средах необходимо использовать различные численные методы. Из-за колебательного характера волновых явлений численные аппроксимации волновых уравнений являются сложными. Аппроксимация волновых задач с допустимой точностью требует использования очень мелкой пространственной дискретизации. Для стандартных численных методов, таких как метод конечных разностей и метод конечных элементов низкого порядка, необходимо сеточно разрешать длину волны посредством десяти ячеек расчетной сетки [3, 4]. В частности, при моделировании задач с высокими частотами требования, налагаемые на расчетные сетки, приводят к задачам с очень большой размерностью, т. е. с высокой вычислительной сложностью.

Классический метод конечных элементов позволяет использовать неструктурированные сетки и проводить расчет в сложных геометрических областях. Данная особенность делает метод чрезвычайно мощным инструментом в различных областях физики и техники при решении различных прикладных задач. Основная трудность классического метода Галеркина [5] при использовании для решения волновых уравнений заключается в том, что аппроксимация уравнения по времени обычно происходит с использованием явных разностных аппроксимаций, что приводит к необходимости инвертирования матрицы масс. Известно, что в классическом методе Галеркина матрица масс не обладает блочной структурой, что приводит к плотной обратной матрице. Данная особенность метода существенно усложняет процесс эффективного численного решения. Предпочтительным для аппроксимации методом для волновых уравнений является разрывный метод Галеркина [6-9]. Данный метод позволяет построить блочно-диагональ-ную матрицу масс, которая может быть легко инвертирована и позволяет построить эффективную вычислительную реализацию. Моделирование распространения волн в упругой среде с использованием разрывного метода Галеркина было исследовано в работах [10-15].

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

Другой особенностью рассматриваемых задач является наличие в грунтах неоднородностей [16, 17]. При этом коэффициенты уравнения могут быть неоднородными, но также в моделируемой среде могут содержаться разномасштабные трещины, которые

приводят к разрыву решений при их аппроксимации с использованием метода LSM (Linear Slip Model). Данное условие на трещине (интерфейсное условие) естественным образом аппроксимируется с использованием разрывного метода Галеркина и учитывается при построении вариационной постановки задачи. Наличие неоднородностей и трещин в грунте существенно сказывается на волновых процессах, и их необходимо учитывать при построении математической модели.

В данной работе мы рассмотрим постановку задачи распространения волн в неоднородной среде. А также представим аппроксимацию математической модели с использованием разрывного метода Галеркина и рассмотрим поглощающие граничные условия. Вычислительная реализация построена с использованием библиотеки FEniCS. Для построения геометрической области используется программа gmsh. Визуализация полученных результатов происходит с использованием программы Paraview [18]. Численные результаты будут представлены для трех модельных задач с неоднород-ностями и трещинами в грунтах.

Постановка задачи

Рассмотрим волновое уравнение, описывающее распределение упругих волн в расчетной области Q:

д2 u

р~дрт - diva (u) = f, t e[0,T], х efi, (1)

где f = f (x,- заданный источник, p - плотность среды, u = u(x,- перемещение в момент времени t.

Уравнение (1) дополняется соотношением между тензором напряжений О и тензором деформаций S

s(u ) = 1 (u + (Vu ) ), (2)

а

(и) = 2ре (и) + XdvvuE, (3)

где Е - единичный тензор, X и ¡л являются параметрами Ламе.

Подставляя соотношения (2), (3) в (1) получим следующее гиперболическое уравнение

Р~~и = рАи + (A + р)graddiv и + f. (4)

Уравнение (4) дополним однородными граничными условиями Неймана

ди (х,t) дп

со следующими начальными условиями

= a, t е [0,T], х еП,

u (x,0) = u0 (x), —( , ) = u1 (x),

, .. п. dt

Указанная выше задача рассматривается в трещиноватой среде Q (рис. 1). Для численного моделирования уравнения упругих волн в трещиноватой среде мы применяем LSM (Linear Slip Model). В частности предполагаем, что толщина трещины пренебрежимо мала. Следуя данной модели для задания трещин, мы имеем линейную зависимость между напряжениями и величиной разрыва в поле перемещения следующим образом

[u ] = Zt, (5)

Рис. 1. Неоднородная область с трещинами

где [и] - скачок поля перемещения на трещине, %=&• п , Z - матрица свойства трещины.

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

Z =

0

0

где 2п и - нормальные и тангенциальные компоненты матрицы свойств соответственно. Конечно-элементная аппроксимация разрывным методом Галеркина

Для аппроксимации задачи по пространственным переменным мы используем разрывный метод Галеркина. Пусть Е - ребро (грань) между элементами К и К2, тогда среднее и скачок вектора и на ребре Е зададим следующей формулой:

и\ = ■

и +и

2

, [и ] =

Пусть Тh триангуляция расчетной области Q и Гh есть множество всех внутренних граней между элементами Тh. Пусть Гc ^ Г h - подмножество всех граней, где поле перемещения является непрерывной функцией и Гf ^ Гh - подмножество граней, представляющих трещины.

Вариационная формулировка волнового уравнения с использованием разрывного метода Галеркина (IPDG, Interior Penalty Discontinuous Galerkin) в трещиноватой среде определяется следующим образом: найти u еVh, такую что

KeThK

я2 д u

2 dt2

vdx +

У К ),е())) - У f{(u )}И ds - У f{r(v )}[u ]

KeThK

ЕеГ ce

Е^ГсЕ

ds

z

z

n

h >

+ у |(А + 2,и)[м][У]ds + ИИds = У \fvdx, V"

ЕеГсЕ EeГfE КеТАК

где у является параметром штрафа и т(и~) = &(и) п .

Для аппроксимации по времени мы используем явную разностную схему:

^ -^^^ + и-+ ^ \ )))- ^ Лт{и")}И^

КеТ ЬК А КеТ ЬК ЕеГсЕ

)

ds +

Yf

EeT CE

+ 2/л) un [v]ds + ^ ^ jz 1 un [v]ds = ^ Jfvdx,

EeT cE

EeT fE

KeT hK

где t - шаг по времени. Поскольку данная схема условно устойчива, следовательно, шаг по времени необходимо ограничивать.

Полученную дискретную задачу запишем в матричной форме

un+1 - 2un + un-1 n Mh-^-+ Khu" = Fh, (1)

где Mh - матрица масс и Kfo - матрица жесткости.

Таким образом, на каждом временном слое мы решаем следующую линейную систему уравнений

Mhun+1 =(2Mh -At2Kh )un -Mhun- + At2Fh,

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

Поглощающее граничное условие

При рассмотрении волновых процессов в ограниченной области необходимо задать поглощающие граничные условия на искусственных границах, чтобы избежать неестественных отражений волн [1, 4, 18, 19]. При проведении расчетов мы будем использовать поглощающее граничное условие первого порядка

ut = -Lan, x (eôQ., t > 0, (8)

где

L =

1/V2Û+I 0 0

0 1/^ 0

0 0 1/^ Граничное условие (8) запишем следующим образом:

1

Vp

-an = Llut, x gDQ, t > 0,

(9)

где

L1 =V

Р

+ A 0 0

0 yfp 0

0 0

Данное поглощающее условие будем налагать в вариационной постановке задачи

и"+1 - 2ип +г

-г-уОх +

д+2

КеТйК ^ ееТье

У>

un+1 - un-1

+

--vdx + У \l 1

At2 ¿¿A 2 At

У\((ип ),e(v ))dx -У\{т(ип )}[v ]d

EeTcE

vds +

KeT hK

Х/Му } ] * + У/ + М\ип + ZÍZ ^ [и ][у № =

ЕеГ СЕ ЕеГ се ЕеГ ге КеТ ЬК

где Гь есть подмножество граней на границе.

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

Mhun+1 +((2Kh -2Mh)u" + Mhun-1 =At2Fh,

а при использовании поглощающих граничных условий имеем

M К+1 + ((2Kh - 2Mh у +

где MdQ,h = J Lluvds.

At

Mh - — M^,h

(10)

-1 = At2 Fh, (11)

0Q

Численные результаты

Приведем результаты численного решения распределения волн в упругой неоднородной среде. Рассмотрим четыре случая:

Задача 1. Однородные коэффициенты.

Задача 2. Слоистые свойства грунта (рис. 2).

Задача 3. Система трещин в однородном грунте (рис. 3).

Задача 4. Слоистый неоднородный грунт с учетом наличия системы трещин (рис. 4).

Распространение упругих волн рассмотрим в области Q = 500 х 500 метров. Для численного решения мы строим расчетную сетку с шагом по сетке h = 4.5 метров, для задачи 1, которая содержит 30016 ячеек и 15007 узлов, для задачи 2 - 33812 ячеек и 17131 узлов, для задачи 3 - 28702 ячеек и 14576 узлов и для задачи 4 - 31092 ячеек и 15635 узлов. Для построения двумерной геометрической модели и генерации сетки воспользуемся программой Gmsh [21]. Для последующей визуализации результатов значения записываем в формате vtk, которые можно визуализировать с использованием программы Paraview [20].

Далее правая часть уравнения дается в виде

f (x, t ) = G (x )P (в)Я (t), где P(0) = (cos(0),sin(0)) - угол источника (0 = 0), R(t) = (l-2a2)e~- волновой

импульс Рикера с a = п f0 52 -

t--

fo

, f0 = 40 Гц, а пространственная функция G (x)

o J

от_I_

Рис. 2. Расчетная область и треугольная расчетная сетка для задачи со слоями. Задача 2

/ У/ / / / * / /

х / / / / ✓ 10.0,01 /V/ / / /

Рис. 3. Расчетная область и треугольная расчетная сетка для задачи с трещинами. Задача 3

/ ' /

У/ . / / /

' УУ / .........../•......../

у /........ /./.

/ г / / / у / у У /

Рис. 4. Расчетная область и треугольная расчетная сетка для задачи с трещинами и неоднородными свойствами. Задача 4

определяется как точечный источник G(х)-5(х-х0), где х0 = (250,250) - центр расчетной области.

Параметры Ламе находятся по следующим формулам:

уЕ Е

% = ■

(1 +у)(1 - 2у),И 2 (1 + у)' где Е = 10^109 - модуль Юнга, V = 0.3 коэффициент Пуассона. В вычислительной области мы

устанавливаем а = 4128.4 [м/c] и в = 1215.3 [м/с], где а —,

скорость p-волны,

в = — - скорость s-волны, р - плотность, равная 2600 [кг/м2]. Для податливости трещин

VP

используем zn = 10 • 10-9[м/Па] и zt = 30 -10"9 [м/Па]. Начальные условия равны нулю и мы вычисляем решение задачи в момент времени T = 0.7 сек с числом временных шагов T = 2900. Величина шага по времени сильно влияет на решение, поэтому расчет

count i

T\max\ 0.7

производится с минимальным шагом по времени т = —-—— =-

Т л 2900

0.00024.

На рисунках 5-8 представлено распределение поля перемещений на различные моменты времени t = 0,15; 0,45 и 0,7. По полученным результатам можно заметить, что в однородной среде волна распространяется равномерно (рис. 5). В случае, когда рассматриваемая область состоит из трех различных подобластей (рис. 6) с разными физическими свойствами, т. е. Е = 10 • 109 Па (верхний слой), Е2 = 5 * 109 Па (средний слой), Е3 = 1 • 109 Па (нижний слой), видно, что волна в верхнем слое проходит через эту область, а в нижнем затухает, так как Е1 > Е2 > Е3 . Как видно из рис. 7, трещины сильно влияют

Рис. 5. Распределение перемещений для задачи с однородными коэффициентами при t = 0.15, t = 0.45, t = 0.7 сек. Задача 1

Рис. 6. Распределение перемещений для задачи с неоднородностью в виде слоев при t = 0.15, t = 0.45, t = 0.7 сек. Задача 2

Рис. 7. Распределение перемещений для задачи с неоднородностью в виде трещин при t = 0.15, t = 0.45, t = 0.7 сек. Задача 3

Рис. 8. Распределение перемещений для задачи с неоднородностью в виде трещин и слоев при t = 0.15, t = 0.45, t = 0.7 сек. Задача 4

на распространение волны, так как они отражают волны. На рис. 8 волны

отражаются от трещин, а также в верхнем слое волна проходит через эту

область, а в нижнем слое затухает, так как рассматриваемая область состоит

из трех различных подобластей с разными физическими свойствами.

Заключение

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

Л и т е р а т у р а

1. CK Bates, DK Phillips, К Grimm, and H Lynn, The seismic evaluation of a naturally fractured tight gas sand reservoir in the wind river basin, Wyoming, Petroleum Geoscience, 7(1):35 44, 2001.

2. GC Cohen and GC Gaunaurd, Higher-order numerical methods for transient wave equations, scientific computation. Applied Mechanics Reviews, 55:B85, 2002.

3. Tatiana Chichinina, Vladimir Sabinin, and Gerardo Ronquillo-Jarillo, Qvoa analysis: P- wave attenuation anisotropy for fracture characterization. Geophysics, 71(3):C37-C48, 2006.

4. Chaur-Jian Hsu and Michael Schoenberg, Elastic waves through a simulated fractured medium.

Geophysics, 58(7):964-977, 1993.

5. Beatrice Riviere, Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation. SIAM, 2008.

6. Wm H Reed and TR Hill, Triangular mesh methods for the neutron transport equation, Los Alamos Report LA-UR-7'3-f79, 1973.

7. Bernardo Cockburn, Fengyan Li, and Chi-Wang Shu, Locally divergence-free discontinuous galerkin methods for the Maxwell equations. Journal of Computational Physics, 194(2):588 610, 2004.

8. Gary Cohen, Xavier Ferrieres, and Sebastien Pernet, A spatial high-order hexahedral discontinuous galerkin method to solve Maxwell's equations in time domain. Journal of Computational Physics, 217(2) :340-363, 2006.

9. Jonas D, De Basabe, Mrinal K, Sen, and Mary F, Wheeler, Elastic wave propagation in fractured media using the discontinuous Galerkin method, GEOPHYSICS, 81(4):T163- T174, 2016.

10. Timo Lahivaara, Discontinuous Galerkin method for time-domain wave problems. University of Eastern Finland, 2010.

11. Yongchae Cho, Maria Vasilyeva, Yalchin Efendiev, and Richard Gibson, Simulation of elastic wave propagation in fractured media with multiscale finite elements. In SEG Technical Program Expanded Abstracts 2016, pages 4003-4007, Society of Exploration Geophysicists, 2016.

12. Marcus J Grote, Anna Sehneebeli, and Dominik Sehotzau, Discontinuous galerkin finite element method for the wave equation, SIAM Journal on Numerical Analysis, 44(6):2408-2431, 2006.

13. Martin Kaser and Michael Dumbser, An arbitrary high-order discontinuous galerkin method for elastic waves on unstructured meshes—i, the two-dimensional isotropic case with external source terms. Geophysical Journal International, 166(2):855-877, 2006.

14. Michael Dumbser, Martin Kaser, and Eleuterio F Toro, An arbitrary high-order discontinuous galerkin method for elastic waves on unstructured meshes-v, local time stepping and p-adaptivity. Geophysical Journal International, 171(2):695-717, 2007.

15. Josep de la Puente, Martin Kaser, Michael Dumbser, and Heiner Igel, An arbitrary high-order discontinuous galerkin method for elastic waves on unstructured meshes-iv, anisotropy. Geophysical Journal International, 169(3) :1210 1228, 2007.

16. Eric T Chung, Yalchin Efendiev, Guanglian Li, and Maria Vasilyeva, Generalized multiscale finite element methods for problems in perforated heterogeneous domains. Applicable Analysis, 95(10):2254-2279, 2016.

17. Eric T Chung, Yalchin Efendiev, Wing Tat Leung, Maria Vasilyeva, and Yating Wang, Online adaptive local multiscale model reduction for heterogeneous problems in perforated domains. Applicable Analysis, pages 1-30, 2016.

18. Utkarsh Avaehit, The paraview guide: a parallel visualization application. 2015.

19. Bjorn Engquist and Andrew Majda, Absorbing boundary conditions for numerical simulation of waves. Proceedings of the National Academy of Sciences, 74(5):1765-1766, 1977.

20. Anders Logg, Kent-Andre Mania!, and Garth Wells. Automated solution of differential equations by the finite element method: The FEniCS book, volume 84. Springer Science & Business Media, 2012.

21. Christophe Geuzaine and Jean-Fran§ois Remade, Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities. International Journal for Numerical Methods in Engineering, 79(11): 1309 1331, 2009.