CC BY
19
УДК 519.6
DOI 10.25205/1818-7900-2020-18-4-66-85
Моделирование 3D волновых полей в неоднородной области со сложной топографией при помощи схемы Лебедева
П. А. Титов
Институт вычислительной математики и математической геофизики СО РАН
Новосибирск, Россия
Аннотация
Численное моделирование широко используется при изучении волновых полей в различных средах. Одним из способов моделирования является разбиение интересуемой области на элементарные объемы и построение конечно-разностной схемы для численной реализации. В работе предполагается, что среда может обладать существенной кривизной поверхности, поэтому задействуется технология построения сетки из криволинейных кубов. Такая сетка обеспечивает хорошую согласованность дискретной и физической моделей среды. Предложен параллельный алгоритм численного решения 3D линейной системы теории упругости, выраженной в скоростях смещений, с использованием криволинейной сетки и явной разностной схемы, созданной на основе схемы Лебедева. Представлены результаты моделирования для неоднородной области. Расчеты проводились с использованием ресурсов ССКЦ СО РАН. Ключевые слова
теория упругости, упругие волны, криволинейная поверхность, криволинейная сетка, суперкомпьютер, моделирование, схема Лебедева Благодарности
Работа выполнена в рамках государственного задания ИВМиМГ СО РАН, проект 0315-2019-0009 (пункты 1, 5-8), а также при поддержке гранта РФФИ № 20-01-00231 (пункты 2-4).
Автор выражает благодарность ССКЦ СО РАН за предоставленные вычислительные ресурсы для проведения численных экспериментов. Для цитирования
Титов П. А. Моделирование 3Б волновых полей в неоднородной области со сложной топографией при помощи схемы Лебедева // Вестник НГУ. Серия: Информационные технологии. 2020. Т. 18, № 4. С. 66-85. БСТ 10.25205/1818-7900-2020-18-4-66-85
Simulation of 3D Wave Fields in Inhomogeneous Domain with Complex Topography Using the Lebedev Scheme
P. A. Titov
Institute of Computational Mathematics and Mathematical Geophysics SB RAS Novosibirsk, Russian Federation
Annotation
Numerical simulation is widely used in the study of wave fields in various media. One of the methods is to divide the domain of interest into elementary volumes and build a finite-difference scheme for numerical implementation. The work assumes that the domain can have a significant curvature of the surface, therefore, the technology of generating a mesh of curved cubes is used. This mesh provides good consistency between the discrete and physical models of the domain. A parallel algorithm is proposed for the numerical solution of a 3D linear system of elasticity theory, expressed via displacement velocities and stresses, using a curvilinear mesh and an explicit difference scheme based
© П. А. Титов, 2020
on the Lebedev scheme. The simulation results are presented. The calculations were carried out using the resources of the SSCC SB RAS. Keywords
theory of elasticity, elastic waves, curved surface, curved mesh, supercomputer, modeling, Lebedev's scheme Acknowledgements
The work was carried out within the framework of the budget project of the ICMMG SB RAS 0315-2019-0009 (sections 1, 5-8), as well as with the support of the RFBR grant No. 20-01-00231 (sections 2-4).
The author expresses his gratitude to the SSCC SB RAS for the provided computing resources for carrying out numerical experiments. For citation
Titov P. A. Simulation of 3D Wave Fields in Inhomogeneous Domain with Complex Topography Using the Lebedev Scheme. Vestnik NSU. Series: Information Technologies, 2020, vol. 18, no. 4, p. 66-85. (in Russ.) DOI 10.25205/ 1818-7900-2020-18-4-66-85
1. Введение
Для лучшего понимания предмета исследований представляем краткое описание процесса: в начальный момент времени полагаем, что в среде не происходит возмущений. В среду помещается источник (заглубленный или поверхностный), который генерирует возмущения. Частицы среды, примыкающие к источнику, начинают колебаться (смещаться) в соответствии с сигналом источника, передавая это возмущение на соседние частицы. В результате этого процесса внутри среды появляется волновое поле. Свойства волнового поля определяются строением и параметрами среды ц - коэффициенты Ламе, р - плотность). Вектор смещения для каждой частицы можно описать через три компоненты в трехмерном случае. Таким образом, восстановить волновое поле - это восстановить вектор смещения для каждой частицы среды в каждый момент времени. На практике также для исследования волновых процессов используется вектор скоростей смещений, зачастую работают именно с ним. В нашем случае численное моделирование волнового поля (векторов скоростей смещений) происходит по дискретному набору точек среды в дискретные моменты времени. Для численного решения задачи задействуется построение конечно-разностной схемы на основе схемы Лебедева, с использованием метода баланса. Отличительной особенностью данной работы является построение криволинейной сетки, хорошо согласованной с формой поверхности области исследования. С некоторыми способами построения криволинейных сеток и их применением для решения актуальных задач можно ознакомиться в [1-3]. В области моделирования волновых полей криволинейные сетки также нашли применение [4; 5], включая работы автора [6-9].
Ввиду размерности задачи и большой вычислительной нагрузки целесообразно использование параллельных технологий для проведения расчетов на многоядерных системах.
В последующем данная работа будет основой при создании 3D параллельного алгоритма моделирования волновых полей в неоднородных средах, характерных для Сибирского региона (например, Баженовская свита).
2. Построение криволинейной сетки
С некоторыми способами построения криволинейной сетки можно ознакомиться в работах [1-3]. Дадим идейное описание процесса построения криволинейной сетки: область, в которой необходимо провести моделирование физических процессов («физическая» область), может обладать значительной кривизной поверхности. Построение криволинейной сетки в этой области делается методом отображений, а именно: между физической областью в пространстве (х1, х2, х3) и областью простой геометрической формы (например, параллелепипед) в пространстве (41, Ч2> Чз) устанавливается взаимно-однозначное соответствие х(ч): Qз ^ Хз. Область в пространстве (Ч1>Ч2>Ч3) называется расчетной. Поскольку рас-
четная область имеет простую форму, ее можно покрыть равномерной регулярной сеткой с ячейками в форме параллелепипедов. В таком случае криволинейная сетка физической области есть результат отображения равномерной сетки в расчетной области. В качестве метода отображения областей был выбран метод трансфинитной интерполяции, как наиболее подходящий для дальнейшей параллельной реализации алгоритма генерации сетки. Дополнением автора является подбор управляющих метрик - функций, позволяющих влиять на параметры отображения и, как следствие, управлять параметрами криволинейной сетки.
В данной работе отличительной особенностью криволинейной сетки является ее ортогональность к свободной поверхности, что накладывает условия на форму поверхностей, с которыми можно работать: если поверхность рассматривать как функцию X 3 = ¿¿(х- > х2) , она должна быть гладкой функцией первого порядка по каждой переменной на участке, на котором она построена для данной области. Пример такой области приведен на рис. 1.
Физическая область Расчетная область
Рис. 1. Пример 3D физической области с топографией поверхности и соответствующей расчетной области Fig. 1. Example of 3D physical domain with free surface topography and corresponding computational domain
Использование такого взаимно-однозначного отображения дает возможность поставленную в физической области задачу в переменных переписать в переменных ( Ч 1 >Ч 2 >Ч з) . Получим видоизмененную задачу, которую требуется решить в расчетной области простой геометрической формы. В расчетной области на регулярной сетке можно использовать известные методы численного решения задачи, например метод конечных разностей.
Опишем формулы для установления взаимно-однозначного соответствия между областями Х(ч)■ (3 — X3 в случае, когда в «физической» области 4 боковые грани и дно области -
плоские и параллельны двум из трех координатных направлений (рис. 1, слева). Для удобства положим, что расчетная область является единичным кубом. Здесь и далее используем обозначения х = (х1 >х2 >х 3) , ч = ( Ч 1 >Ч 2 >Ч 3 ) .
Отображение граней областей считаем заданным. Оно также строится
методом трансфинитной интерполяции, только с применением двумерных формул. Подробней можно ознакомиться с этими формулами в работе [3].
Имея отображения граней областей друг на друга, можно эти формулы интерполировать внутрь, получив тем самым полное отображение областей. Для управления свойствами сетки используются переходные функции:
«1, о ( Ч1) = ( 1 + 2 Ч1) ( Ч1 - 1) 2, а£0( Ч1) = 1 - < О (Ч1 ) , а1,1 (Ч1) = Ч1 ( 1 - Ч1 ) 2, 4 1 ( ЧО = ( Ч1 - 1) Ч12,
а2, о ( Ч2) = ( 1 + 2 Ч2 ) ( Ч2 - 1) 2, а|о(Ч2) = 1 - а£ о(Ч2 ) , а£ 1 (42) = Ч2 ( 1 - Ч2) 2, а2, 1 ( Ч2) = ( Ч2 - 1)Ч22,
«I3 о ( Чз) = ( 1 + 2 Чз ) ( Чз - 1) 2, а з,о(Чз ) = 1 - а?,о ( Чз ) , а£ 1 ( Чз ) = Чз ( 1 - Чз ) 2*,
а3, 1 ( Чз) = ( Чз - 1) Чз2 * (1)
Здесь параметр 2k - четная степень, является нововведением автора. Этот параметр подбирается отдельно для каждой физической области. Размер области и кривизна поверхности будут играть основную роль при выборе значения k. Также считаем, что криволинейная поверхность является гладкой функцией второго порядка . Для проверки
перехлеста ячеек используется Якобиан / = й е t (-—( х ¿) ) . Во всей области должно выполов у
няться условие . Здесь и далее используем обозначения для производных
¿(Л, ¿ = и, 3.
Теперь более подробно опишем итоговые формулы для установления взаимно-однозначного отображения между единичным кубом (расчетной областью) и физической областью
Х(Ч)^Х з:
Ч2. Чз) = Ц2. Чз) +
и1,о(.Чз)ШЧ1.Ч2.Ъ) ~ $2ХЧ1.Ч2.®)~\ + <4,\.(Яг)[хща(Я1,Чг>0) ~ ^+
Я!,о(<7З)1Ж<71.<72. 1) -S2i.q1.q2.1)] + <4,\.(Яг)[х1Ца(Я\,Чг. 1) ~ ^ Ч 1 , Ч 2,1 ) ] , ¿ = 1, 2, 3, (2)
где используются обозначения
^Ч Ч 1 ,Ч 2 ,Ч з ) = а 1 о ( Ч 1 ) х ¿( 0 ,ч 2 ,Ч з) + а £ 1 ( Ч 1 ) х ^ ч Д 0,ч 2 ,Ч з) + , «2,0(^1)^(1.(72,(7з) + с4д((71)хг^(1,С72,/ = 1,2,3,
^ <7з) = <7з) +
О, (7з) -5/^(67!, 0,(73)] +
«2,о(<72)1Ж<71. Мз)] 1. Чз) ~
Ч 1 ,1,Ч з) ] , I = 1, 2 , 3 .
3. Постановка задачи 3.1. Постановка задачи в «физической» области
Процесс распространения волновых полей моделируется при помощи численного решения полной линейной 3D системы теории упругости, записанной в скоростях смещений и напряжениях. Физическая область в общем случае обладает криволинейной свободной поверхностью (при этом, как в параграфе 2 данной статьи, считаем, что поверхность есть функция, удовлетворяющая условиям гладкости), среда - изотропная, не однородная. Граница области разделяется на две части: - свободная криволинейная поверхность и - четыре боковые и рань, противоположная свободной поверхности. В «физической» области система принимает вид
Р Ус ^ = ¿Г + 4 (^12) + Кз)
^ I М = £ (*12) + £ (*22) + ОЪз)
^ ^ ^ ^
^(из)=^(а1з)+^(а2з) + з^(сГзз) ^ ^ ^ ^
- (ап) = (Я + 2ц) — (щ) + я — (и2) + я — (и3)
Т2 2 ) = Я 0 + (Я + 2 ^¿С и 2) + Я^-С и з) (3)
^(азз) = + + (Я +
д д д ^ д д ^ ^ д
Начальные условия при t = 0:
и ¿(Х>0 ) = 0 > £( и ¿) (х>0 ) = 0 (4)
Граничные условия на 3 Г:
и I | аг = 0> I = 1> 2 > 3 (5)
Условия на свободной поверхности 3 5:
"1011 + "2012 + п3а13 = 0, ща12 + П2а22 + п3а23 = 0, пхахъ + п2а23 + п3а33 = 0> V t (6)
Здесь ( и - (х> ^ ) > и2 (х> ^) > из (х> ) ) - вектор скоростей смещения вдоль координатных осей Ох-> ОХ2 > Охз, а = ( (¿у) - симметричный тензор напряжений (оуу = (Ту ¿), р (х) - плотность, (х) - скорость поперечных волн в среде, (X - скорость продольных волн в среде. Я(х)
= ( — 2 1^2) /р , д (X = ^2/Р - параметры Ламе, 35 - свободная поверхность,
3 Г - граница области, не включая свободную поверхность, [п- > П2 > П3 ] т - единичный вектор нормали к свободной поверхности .
3.2. Перенос исходной задачи на «расчетную» область
Сформулирована математическая постановка задачи в «физической» области при помощи уравнений (3)-(6). Используя правила замены переменных, перепишем эти уравнения в терминах . Учитывая замену , получаем
и (4) = и(х- (Ч) >Х2(ЧО > х 3(4) Л )>1 = 1> 2> 3.
Необходимо сохранить дивергентный вид первых трех уравнений для дальнейшего применения метода баланса, описанного в [10]. Аналогично работе [14], система (3) преобразуется в систему вида
jp lt(u2)=^г(а4)+4+iм
a a _ a _ a _
3 3
a a v- a
¿=1 ¿=1
Z3 a a
(7)
3 3
a
L3 dqt
i=1 t=l
Z3 a a
¿=1
3 3
a a v
¿=1 ¿=1
3 3
v- a a v-
¿=1 ¿=1
3 3
a v- a a v-
i=1 t=l Начальные условия при t = 0:
77,Гг m = П
dt
Граничные условия на д Г:
a a
dx1 d4i
a a
dx1 dqt
a a
dx1
ui(x, 0) = 0, £( и¿) (х, 0) = 0 (8)
u i | аг = 0, i = 1, 2, 3 (9)
не меняются.
На свободной поверхности д5 условия видоизменяются:
cfp=0, оГ=0, ôçT=0 (10)
Здесь
=/2 3= 1 <1)1 , ^2 = /!f= 1 <2)1 , ^3 = /!f= 1 Чз)Ъ ь
3 а 3 а 3 а
i=l i=l i=l
Оп
3 3 3
t=l t=l t=l
t—lr ^ ^ м ^ г—lr ^ ч ^
Ч1) =/" 1 [а^(х2)а^(хз)—а^(хз^^ ) ] , Ч2 ) =/- 1 ^ )^(хз ) — ] ,
в ¿( Чз) = /" 1 [¿(х2)4(хз) —¿(хз )£ (х2) ] , ¿( Ч1 ) = /- 1 ¿(хз )а4(х 1 ) —
¿(х1 )4(хз) ] ,
в ¿( Ч2) = /- 1 [¿(хз)^(х!) — ¿(х1)^(хз) ] , ¿(Чз ) = /- 1 [¿( х з)^(х 1 ) —
¿(х1 )^(хз) ] ,
в ¿( Ч1) = /- 1 [¿(х 1 )^(х2 ) — ¿(х2)^(хз) ] , ¿(42 ) = /- 1 [¿(х1)^(х2 ) — ¿^¿(х1 ) ] ,
-а-( ч з) = / - 1 [-а-(х1 ) -— (х2) — —— (х2 )та-(х1) ] , д
] = <1е1{— (хО) дц]-
д д д д д д
д д д дЦз о <71 дц2 х1) ( х2) — ( х з) — — ( х 1) — (х2 ) — (х з ) — — ( х 1) — ( х2 ) — (х з ).
Далее для системы (7) с граничными и начальными условиям (8)-(10) опишем созданную конечно-разностную схему.
4. Численная реализация с использованием явной разностной схемы Лебедева
В данном параграфе предлагается явная по пространству и времени конечно-разностная схема для численного решения задачи (7)-(10) внутри «расчетной» области. Как и в работе [11], задействуем принципы расположения элементов в ячейке аналогично схеме Лебедева. На рис. 2 изображена стандартная кубическая ячейка ( / >у > к ) в расчетной области с расположением в ней информации. В узлах, обозначенных треугольными маркерами, хранятся компоненты вектора скоростей смещения (и 1>и 2 >и з ) , в узлах, обозначенных круглыми маркерами, хранятся компоненты тензора напряжений . В целях удобства индексации при реализации разностной схемы, все узлы, содержащие скорости смещений, были разделены на группы , аналогично для напряжений разделение происходит на группы .
Для получения большей точности при расчетах все метрические коэффициенты и якобиан рассчитываются на сетке с половинным шагом. Как следствие, для них предусмотрена «удвоенная» индексация. Например, узлу будут соответствовать
/2 I, 2/, 2 — ( х^)) 2 I,2 у,2 к узлу 4 ( ¿>)> к ) - /2 1 - 1, 2/ - 1, 2 к - 1 > — (х:5:)) 2 I- 1, 2у - 1, 2 к - 1 >5>ГП = "Чт "Чт
1, 2 > 3 .
Параметры среды р^ , у, Я ^ , у, ц¿,у, к располагаются в центрах ячеек, т. е. в узлах 4 '( / ,_/', /) .
Для удобства считаем пространственный шаг сетки одинаковым по всем направлениям Л = Л/ ± = Л/ 2 = Л/ 3 (для метрических коэффициентов и Якобиана пространственный шаг берем Л/ 2 ), временной шаг обозначим т. Вся расчетная область покрывается регулярной квадратной сеткой М X Ь X К ячеек ( 2 М + 1 X 2 Ь + 1 X 2 К + 1 для метрических коэффициентов и Якобиана). Сетка узлов для каждой группы имеет размерность: { 1 } — (М + 1 X Ь X К) , { 2 } — (М X Ь + 1 X К) , { 3 } — (М X Ь X К + 1 ) , { 4} — (М X Ь X К) , { 1 '} — (М X Ь + 1 X К + 1 ) , { 2 '} — (М + 1 X Ь X К + 1 ) , { 3 '} — (М + 1 X Ь + 1 X К) , {4'} — (М + 1 X Ь + 1 X К + 1 ) .
Для произвольной функции А ( б/, Г,*) обозначим А [у к = А ( I Л,/ Л, кЛ, пт) , а для Якобиана и метрических коэффициентов . Введем обозначения разностных аналогов дифференциальных операторов:
О* [А ] Пм/2 = т "1 (А52 — А »,/) = А * ( ¿Л,уЛ,к Л, (п + 1 / 2 ) т) ,
О * [А ] [у, / = т- 1 (А П5/2—А П;,1/2) = А * ( ¿Л,у Л, к Л, пт) ,
Я [А] П- 1/2,;,к = Л- 1 (А [у к — АП- 1 ,у к) = А ^ (¿Л — Л/2,уЛ, кЛ, пт) ,
0 2 [А] [у-1 / 2, к = Л- 1 (А [у к — АПу- 1, к) = А ^(¿Л,уЛ — Л/ 2 , к Л, пт) ,
= = АЧг{1к,]к,кк-к/2,пх)
01 ] 2 ¿, 2у 2 к = Л- 1 2 ¿+1, 2у 2 к — Я2 ¿ - 1, 2у 2 ^ = ЯЧ1(¿Л,У Л, к Л, пт) ,
О 2 ] 2 ¿, 2у 2 к = Л- 1 2 ¿, 2/+1, 2 к —Я 2 ¿, 2/ - 1, 2 к) = Я <?2( ¿Л,У Л, к Л, пт) ,
О 3 ] 2 ¿, 2у 2 к = Л- 1 2 ¿, 2у 2 к+1 — Я 2 ¿, 2у 2 к - 1) = Я <?3 ( ¿Л,У Л, к Л,пт).
Рмс. 2. Расположение элементов в ячейке (г, j, k) Fig. 2. Elements' positioning in (г, j, k) cell
Далее вводим обозначения элементов для разных групп: - элементы тен-
зора напряжений группы (1'), I>ж = 1, 2 > 3 и соответственно для <тр ( 1 ') =
= з= 1 "д^Х4 1) <1 1( 1 '), как и для остальных О^ — Од" ; и ( 1 ) ¿у, / - элемент вектора
скоростей смещений группы {1}. Аналогично обозначаются элементы групп {2'}, {3'}, {4'} и {2}, {3}, {4}.
Тогда разностная схема для (7) примет следующий вид. Для скоростей смещений
= /211-1,2у,2*{01 № (.4%-ïllj,k+D2 [ат (3')]-;+11//2,А
+ D3 [dT(2')
= Î2i-l,2j,2k{^ + D3 W (2')
= J2i-l,2j,2k{P
+ D3 [dT(2')
= hu2j-l,2k{P
+ D3 [dT(l') Pi,y-l/2,fcôt[W2(2)]^Î/2
= J2i12j-l,2k{D + D3 W (1')
= J2i12j-l,2k{D + D3 № (1')
Pi,j,k-i/2Dt[Ui(.3)]tH/2
= h^jM-iiD + D3 [dT(4')
= J2i,2j,2k-l{D + D3 [dT (4')
n+1/2 i,y,fc+l/2J
n+1/2 J
тП+1/2
i,j,k+l/2J
[*r mTlILk+^2 [dr
n+1/2 -, t,y,fc+l/2J
qn+l/2
n+1/2 -,
i,y,fc+l/2J
лп+1/2
[<T (3')]^y,fc + D2 [dT (4')]i'/-i/2,i
n+1/2 -, i,y,fc+l/2J
-.n+1/2
n+1/2 -, i,y,fc+l/2J
тП+1/2
n+1/2 -,
i,y,fc-l/2J
n+1/2
n+1/2 -, t,y,fc-l/2J
qn+l/2
Ри,к-1/20ЛщШ1Ц/2
= Jlàm-i{Di № + D2 [<*Г (!')]■
+ D3 [aT
= /2l1-l,2y-l,2fc-l(ßl [<*T d')]^ + D2 (2')] -AM, + № Pi- 1/2J-1/2,к-1/2Dt [u2 (4) ] ПцЦ2
= J2i1-l,2j-l,2k-l{Dl W + D2 [C (2%ll-l2/2,k
+ D3 [аГ C3')]^i/2} Pi-l/2,j-l/2,k-l/2Dt[U3(.^]lÜ/2 = J2i-l,2j-l,2k-l{Dl W (1 0]^/2J,fc + 02 [<*T (2')]-/ЛА, + D3 W (3')];:I^21/2} (H)
И для напряжений
I i,j,k
— (Я + 2^)jj_1/2,A;-i/2{i7lx12t,2y-l,2A:-l^l[lil(4)]f+i/2j,A;
i,j—1/2,к "Г"
iD3 [Щ (2)]
i,j,k-1/2) +
^■i,) — l/2,k—l/2{.4l x2 2i,2j—1,2k —1^1 [^2 (4)] ¿+l/2j',fc + ^2X2 2i,2j-l,2k- iD2 [w2(3)] i,j—1/2,к + ^3 Xi 2i,2j-l,2k- i,j,k-1/2} +
1/2,k—l/2{.4l x2 2i,2j—1,2k —1^1 [^3 (4)] ¿+l/2j',fc + ^2X2 2i,2j-l,2k- i,j—1/2,к
+ ^3 Xi 2i,2j-l,2k- iD3[u3(2)] i,j,k-1/2}
= (Я + 2^)j_1/2j',A;-i/2{i7lx12t-l,2y,2A:-li)l[lil(3)]f-i/2j,A; + X! 2i—1,2j,2k— 1^2 [U1 (4)] i,j+l/2,k "T"
^ЗХ! 2i—1,2j,2k —
l№l(l)]
i,j,к—1/2} +
1^2 ["2(4)] i,j+l/2,k
+ 9з Xi 2i—1,2j,2k— i,j,k—1/2} +
Я-i-1/2,У,A:-1/2{^lXz 2i-l,2y,2fc-l£)l[u3(3)]f-i/2,y,A; + <?2 x2 2i-l,2j,2k- iD2[u3m i,j+l/2,k
+ 9з Xi 2i—1,2j,2k— iD3[u3(ï)] i.j.k-lß} Dt[ai i(3')]î!y,fc
= (Я + 2/l)j_1/2,y-l/2,A;{i7lx12t-l,2y-l,2A:i)l[lil(2)]f-i/2,y,A;
i,j—1/2,к "Г"
^ЗХ! 2i—l,2j—l,2k^3 [^1 (4)] ¿^y,fc + l/2} + ^■i-l/2,j-l/2,ki.4l x2 2i —1,2 j—l,2k^l \М2 (2)] ¿^-l/2,y,fc + ¿72 x2 2i-l,2j-l,2ku2 ["2(1)] i,j—1/2,к + i73x12i-l,2y-l,2A:i)3[u2(4)]"y,A; + i/2} +
^г-1/2,у-1/2,/с{'71 х2 21-1,2у-1,2&^1[из(2)]Г-1/2,у,А: + 42 х2 21-1,2 )-1,2к^2 [м3(1)]^у-1/2,А: + Чзх1 21-1,2у-1,2А:^з[из(4)]"у,А;+1/2}
^[^11(40]";,л
= (я + 2д);у 1к{ц1х1 21,2у,2/с^1[м1(1)]г+1/2,у,а:
+ 2г,2У,2/с^2 [и1 (2)]"у+1/2,/с + 21,2У,2/с^з[и1(3)]1,у,А; + 1/2} + ^,],к{Ч1х2 21,2],2к^1 [^2(1)] 1+1/2,]',к + ^2 х2 2£,2У,2/с^2 [иг(2)] 1,]+1/2,к + 9з 21,2У,2/с^з[и2(3)]^у,А:+1/2} + '^¿,У,А:{(71х2 21,2У,2/с^1[из(1)]Г+1/2,У,А: + ^2 х2 21,2],2к^2 [из(2)]^у+1/2,А: + Чзх1 21,2 ],2к^3 [изШ,1к+1/2} (12)
Разностные аналоги для остальных уравнений на компоненты тензора напряжений ( (¿у) для групп определяются аналогичным образом.
Начальные условия:
и0 I , у, к ( 1) = и0 ¿, у, к ( 2 ) = и0 I , у, /( 3 ) = и0 I , у, к ( 4) = 0, 5 = 1, 2, 3 (13)
остаются неизменными.
Граничные условия на :
^(1)1, у, к = и5(1)м+1,у,А: = и5(1)и,А: = = и$(1)^уд = 0'
^(2)1, у, к = из(?-)м,],к = ^(2)1,1, к = из(.2)ц+1к = и5(2){,уд = 0, ^(3)1, у, к = из(^>Ум,},к = к = из(^>У1,1, к = =
"5(4)1,;,Л = и5(4)£+1,м = и*(4)£и = = И5(4)£ул+1 = 0,5 = 1,2,3,п > 0,
^лП+1/2 _ ^^/чП+1/2 _ _ _ р.
) 1,],к ~ )м+1,],к ~ )1,},к ~ )м+1,],к ~ и<
Г1/лп+1/2 _ Г1,лп+1/2 _ гг.,лп+1/2 _ гг.,лп+1/2 _п
)1,1,к ~ )1,1+1,к ~ )1,1,к ~ Л,1+1,/с ~ и<
<Ъп( 1')Ц,к+1 = °1т(2') 1=0- 1т = 1-2,3, П > 0 (14)
Граница свободной поверхности проходит вдоль граней ячеек с индексом . Как
и в работе [9], здесь задействуется «фиктивный» узел при реализации конечно-разностной схемы для групп .
На свободной поверхности условия примут вид
Л (Т) п*[/2 = о, «г сю :;;/2 = о,*г п; = 0,
л-12') ^1/2 = о,ог(2') = о,^(2') ;;;<2 = п.
^(30 ^ = 0,^(3') = С» = а <Г № = 0>°т М = 0>«Г М ££ = 0 >п > о. (щ
В (11) распишем для примера разностный аналог для 01 ( 4 ' ) /2:
( Ч0] 2 ¿, 2 У, 2 3 ( 4 0^/2}. Остальные — в наборах { 1 } — {4 } получаются аналогичным образом. Параметры среды на гранях ячеек определяются формулами:
Р£-1/2,у„к = (Рц,к + Р1-1,у,/с)/2, Рц-1/2,к = (Рц,к + Рц-1,к)/2. Р1,},к-1/2 = (Рц,к + Рц,к-1)/2,
Р1—1/2,} —1/2,к—1/2
= (Р1,],к + Р1-1,],к + Р1,)-1,к + Р1,),к-1 + Р1-1,)-1,к + Р1-1,),к-1 + Р1,)-1,к-1 + Р1-1,у-1,/с-1)/8'
(я + 2д);_1/2,у_1/2,а:
= 4 [(я + 2/0^* + (я + 2/0ГЛ,/,л + (я + 2 /0 + (я + г^гл.у-и]-1. (я + 2д);_1/2,уА;-1/2
= 4 [(я + 2/0 • ^ + (я + 2/0гл,/,л + (я + + (я + г/огд,/,*-!]-1' (я + е,у—1/2,/с—1/2
= 4 [(я + 2/О^л + (я + 2 д)+ (я + 2 д) + (я + 2/0^,^г1, 1/2, у-1/2, к = + + ^1,]-1,к + У-1-1,]-1,к\ Х>
^1-1/2,),к-1/2 = + + ^1,],к-1 + Й-1,У,/с-1]
\11,)-1/2,к-1/2 ~ + ^,]-1,к + ^,]',к-1 + Й,у-1,/с-1]
Разностные формулы для метрических коэффициентов:
/21,2У,2/с = [х1]21,2],2к^2 [-^2]21,2],2к^3 [хз]21,2У,2/с
+ ^2[х1]21,2У,2А:^з[х2]21,2У,2А:^1[хз]21,2У,2А: + 03 [х1]21,2У,2к^1 [хг]21,2],2к^2 [хз]21,2У,2/с — ^3 [х1]21,2],2к^2 [хг]21,2У,2/с^1 [хз]21,2У,2/с — ^2 [х1]21,2],2к^1 [х2]21,2],2к^3 [хз]21,2У,2/с — ^ 1 [х1]21,2],2к^3 [хг]21,2],2к^2 [хз]21,2У,2/с<
[^((7З)]21,2У,2А: = ¡21^2],2к(^1 [х1]г1,2У,2А:^2 [хг]21,2У,2А: — [хг]21,2У,2А:^2 [Х1]21,2У,2А:)-
Для остальных ( Чт) , 5, ш = 1, 2 , 3 , разностные аналоги получаются таким же образом. ох5
Опишем отдельно реализацию источника возмущений в среде. В данной работе используется источник типа "центр давления" = ( 4, £) = (Х1 ( 4) , х2( Ч) ,хз( Ч ) , 0, 5 = 1, 2 , 3.
Источник помещается в центр ячейки, в узел ( /,_/', // ) , принадлежащий набору {4 '}. Для численного моделирования возмущений в среде на каждом временном шаге необходимо обновлять значения переменных, отвечающих за главные напряжения вблизи источника, а именно:
л
1
«u(l')S/2
iwYïif ,
«ъсюгйй =
iF3(4'),;,i;£/2. «u(2')S/2 :
Ï^Jj/2 ,
îWC/2,
;f2(4')")i/2 ,
«bPOÏÏïS =
i^v2.
«иСЮ^+^С^
п,лп+1/2 1 F гл,лп+1/2
n+1/2 1рГД,.п+1/2
n,.n+1/2 lp n+1/2
п,лn+1/2 1 F алП+1/2
<x22(2')^1/2+iF2(4')-;,)f
п,лп+1/2 1 F f.,.n+l/2
022 +4^2(4 j Uj k
n,.n+1/2 1 F гл,лп+1/2 ^ЗЗС^ +4^3(4 )i jk
0iid = 0h(I');:ïï; +
n+1/2
=0n(i')^î +
n+1/2
022(1= 022(1')^ +
n+1/2
ar.n+1/2 __/чn-t-i/z .
Л,У+1,/с+1 _ 022v-^ JjJ+i.fc+i +
n+1/2
«»ffffl = »»ciorSfi +
ai^n+1/2 __/ул-t-i/z .
Л,У+1,/с+1 _ 03311 Jj,y+i,fc+i +
n+1/2
011(20^5=+
n+1/2
(-nhJl+1/2 __7i-t-±/z .
0"lllz Jt + lj,fc+l - Ji+lj,fc+l +
n+1/2
022(20^5=^(2')^;;+
n+1/2
sr.i.n+1/2 __^/yfi-t-i/z .
°221/ Jt + lj,fc+l _ a22lZ Jt+lj',fc+l +
n+1/2
= »»PO&'fi +
t+lj,fc
^лП+1/2 __s^/.n-t- i/z .
0"33lZ jt + l,y,fc+l - °331/ jt+l,y,fc+l +
n+1/2
n+l/2 _
Го/чП+1/2 lp
°ll(3 +4Ы4 )i,j,k
«ni* OïSfë = «nPOïSfô +
i+l,j,k
.n+l/2
". 1 ( 3 ') = 1 (3 ') „Ä . = "i.^'JÎu +
^n+1/2 !,k
= <x22(3')S/2+iF2(4')-f
t+l,y,fc
Nn+1/2
!,k
^о/чП+1/2
22 (3 ),,+;,.
ïf2(4')u
n+l/2 к
rr,n+l/2 1 р г.,лп+1/2
n+l/2
^зз(30Г++й = -зз(З')^; +
n+l/2
i,j,k
\wfCiT .
"3a (30^fi = "33( 3 ') ÏÏÎS +.f3(4')7.7i/2 . "зз(3')?++^и = «за(30^„ + ^3 (4')"li/2 .
n+l/2
^n+1/2 i,i,k
w 'TCiT
n+l/2 n+l/2
^2(4 о-;,
i.i.k
i.i.k
^2 (402 ,
=-33(40^r+^(40-If
n+l/2
1/2 i,i,k
h^itiT
= "пС^Г +
n+l/2
= tfnwnr +
n+l/2
= «aC^ET +
<r33(4')""/2 = ff33(4')"^/2 +
n+l/2
Представленная схема ( 11)—( 15) обладает вторым порядком точности по пространству и времени во всей области, кроме границ раздела двух сред; на них схема обладает первым порядком точности по пространству и вторым порядком точности по времени. Критерием
устойчивости схемы является видоизмененное условие Куранта л/ЗУр тахТ < Jmiпh, где J min - минимальное значение Якобиана во всей области, Vp mах - максимальное значение скорости продольной волны во всей области. Физический смысл условия - за один временной шаг т волна не должна продвинуться по среде дальше, чем диаметр самой маленькой ячейки (максимальное расстояние между двумя точками ячейки).
5. Программная параллельная реализация алгоритма
Для решения данной задачи необходимо применение суперкомпьютерных вычислительных систем и технологий параллельного программирования. Для вышеуказанной конечно-разностной схемы (11)—( 15) был разработан 3D параллельный алгоритм и его параллельная программная реализация при помощи языка Fortran 95, а также с использованием технологии MPI [12]. Также автором отдельно был реализован параллельный алгоритм генерации 3D криволинейной сетки, также на языке Fortran 95 с применением MPI.
Суть параллельного алгоритма заключается в декомпозиции всей расчетной области на отдельные подобласти меньшего размера (кубы), в которых вычисления производятся одновременно. Каждой подобласти назначается вычислительный поток (параллельный процесс). При необходимости производится обмен данными между соседними процессами. Топология связей между процессами задается пользователем при помощи встроенных средств MPI.
В данной работе используется топология связей "3D-куб", внутри которой у каждого процесса есть координаты (рис. 3).
Рис. 3. Топология обменов данными между соседними процессами Fig. 3. Data exchange topology between neighboring processes
Ввиду особенностей конечно разностной схемы между соседними процессами необходимо производить обмен данными на каждом временном шаге (пересылки). Используются неблокирующие команды MPI_ISend, МР1_Жесу для пересылок. Это позволяет оптимизировать время работы программы. Для возможности реализации таких обменов данными, вычисления внутри каждой подобласти были разделены на 2 части: вычисление в приграничных элементах куба и внутри куба. Каждая итерация по времени состоит из последовательных этапов: сначала для каждого процесса производятся расчеты напряжений и скоро-
стей смещений в приграничной зоне, затем вызываются команды обмена данными MPI_ISend, МР1_Жесу, и следом продолжаются расчеты напряжений и скоростей смещений во внутренней части куба.
Для численного моделирования распространения волн рассматривая модель неоднородной среды, срез которой представлен на рис. 4. Срез сделан параллельно плоскости Х2, по оси Y срез содержит точку у = 1. Свободная поверхность выражена функцией:
По форме поверхность схожа с представленной в работе автора [9], с измененными параметрами относительно ширины и высоты криволинейного участка.
«Физическая» область покрывается криволинейной сеткой, локально-ортогональной возле свободной поверхности (см. рис. 4). Именно такой способ построения позволяет удобно и более точно реализовать условия на свободной поверхности (15). Сетка конструируется с использованием формул (2), для управляющих функций (1) выбран параметр k = 10. Модель состоит из 4-х слоев, каждый из которых характеризуется параметрами: слои I и III - р = 1000 к г /м 3 , Кр = 4000 м / с, = 2 000 м / с, слой II - р = 800 к г /м 3 , Кр = 22 3 6 м / с, = 1 1 1 8 м / с, слой IV - р = 900 к г / м 3 , Кр = 3 1 6 2 м / с, = 1290 м / с .
Источник волнового поля - точечный, типа «центр давления». Координаты источника (х, у, z) = ( 0 . 2 5 к м , 1 . 0 к м, 0 . 2 5 к м ) . В роли временной составляющей сигнала выступает импульс Рикера R( t) = ( 1 — 2 7Г2ш>2( t — t0) 2) е_ 27r2ûj2 (t_ {о)2 с несущей частотой ш = 30Гц и временным сдвигом to = 0 . 0 5 секунды. Тогда функции, отвечающие за работу источника, принимают вид: Уо ,z — z o ) R( t) .
Размерность тестовой задачи - 1000 х 1000 х 501 ячеек для переменных и для параметров среды, 2001 х 2001 1003 ячеек для криволинейной сетки и связанных с ней метрических коэффициентов. Время моделирования процесса - 1.4 секунды, временной шаг At =
= 0.0001 секунды.
6. Модель для расчетов
„ , , 3 Лп(х-1\ , 1 Л2п(х-1\,,3 ,4п(у-1\ , 1 -
Z = 1 + (— cos(—-—-) +— cos(--—cos(—-—-) +— cos(-
Чб v 5 J 16 v 5 JJK\b v 5 J 16 v
x < 1.4,0.6 < y < 1.4,
z = 1 иначе.
IV
0.8
06
Рис. 4. Модель неоднородной среды (срез в плоскости XZ) и соответствующая локально-ортогональная сетка
и
Fig. 4. The Model of heterogeneous medium (XZ section) and corresponding locally-orthogonal mesh
1.2 1.4
16
18
x
7. Результаты численного моделирования
Расчеты проводились на кластере НКС-1П ССКЦ СО РАН (www2.sscc.ru). Каждый узел включает в себя два 16-ядерных процессора Intel Xeon E5-2697v4 (2.6 ГГц) и 256 Гб оперативной памяти. Ввиду высокой требовательности задачи относительно оперативной памяти, было задействовано два узла. Общее время вычислений составило около 18 часов. На рис. 5 представлены результаты численного моделирования.
0-2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
Рис. 5. Снимки волнового поля Fig. 5. Wave field snapshots
Как можно наблюдать, отсутствуют дифракционные волны при отражении от поверхности. Однако наблюдаются другие эффекты, связанные с применением схемы Лебедева и криволинейных сеток при реализации конечно-разностной схемы. На снимках отображение амплитуд волнового поля сделано сильно контрастным, чтобы была возможность наблюдать низкоамплитудный дефект, отмеченный на снимках красной стрелкой. Можно наблюдать, как коническая волна опережает фронт Fp-волны. С подобным сталкивались и авторы работы [11], применяя схему Лебедева при решении 3D задачи моделирования волн в анизотропной среде. Однако ввиду небольшой амплитуды артефакта с ним научились работать, и он не доставляет явных неудобств при практических задачах.
8. Заключение
В работе был представлен параллельный алгоритм моделирования 3D волновых полей в изотропных средах, явная конечно-разностная схема для численной реализации и результаты численного моделирования. Новизна работы включает: параллельный 3D алгоритм реше-
ния динамической задачи геофизики с использованием схемы Лебедева в обобщенных координатах, параллельный алгоритм построения 3D криволинейной сетки, новый способ аппроксимации граничных условий на свободной поверхности с участием такой сетки.
К положительным качествам алгоритма можно отнести его способность устранять дифракционные волны при отражении от поверхности, которые характерны при дискретизации области на стандартные прямоугольные ячейки [13]; также его применимость в областях, имеющих жидкостные включения (например, резервуары с нефтью). К отрицательным - высокая требовательность к оперативной памяти и проявление не физичных эффектов в виде фиктивных конических волн. Оба отрицательных фактора обусловлены использованием схемы Лебедева и обобщенных координат (с точки зрения уравнений получается аналог анизотропной задачи). Фиктивные волны легко выделяются и не представляют проблем при анализе результатов моделирования. В дальнейшем планируется реализация поглощающих граничных условий PML [14] и оптимизация алгоритма по объему требуемой оперативной памяти для применения к задачам реального масштаба.
Список литературы
1. Liseikin V. D. Grid generation method. 2nd ed. Berlin, Springer, 2010. ISBN 978-90-4812912-6
2. Vaseva I. A., Liseikin V. D., Likhanova Yu. V., Morokov Yu. N. An elliptic method for construction of adaptive spatial grids. Russ. J. Numer. Anal. Math. Modelling, 2009, vol. 1, p. 65-78. DOI 10.1515/RJNAMM.2009.006
3. Лисейкин В. Д. Разностные сетки. Теория и приложения Новосибирск, 2014. 245 с.
4. Appelo D., Petersson N. A. A stable finite difference method for the elastic wave equation on complex geometries with free surfaces. Communications in computational physics 5, 2007, no. 1, p. 84-107.
5. Jose M. Carcione. The wave equation in generalized coordinates. Geophysics 59, 1994, no. 12, p. 1911-1919. DOI 10.1190/1.1443578
6. Титов П. А. Алгоритм и программа моделирования 2D-волновых полей в областях с криволинейной свободной поверхностью // Материалы конференции "Научный сервис в сети Интернет - 2014". Новороссийск, Абрау-Дюрсо, 2014. C. 446-455.
7. Титов П. А. Моделирование упругих волн в средах со сложной топографией свободной поверхности // Вестник НГУ. Серия: Информационные технологии. 2018. Т. 16, № 4. С.153-166.
8. Titov P. A technology of modeling of elastic waves propagation in media with complex 3D geometry of the surface with the aid of high performance cluster. Russian Supercomputing Days, Proc, 2016, p. 1020-1031.
9. Titov P. The Simulation of 3D Wave Fields in Complex Topography Media. In: Russian Supercomputing Days, CCIS 1129, 2019, p. 451-462. DOI 10.1007/978-3-030-36592-9_37
10. Vishnevsky D., Lisitsa V., Tcheverda V., Reshetova G. Numerical study of the interface errors of finite-difference simulations of seismic waves. Geophysics, 2014, vo. 79 (4), p. 219232. DOI 10.1190/geo2013-0299.1
11. Vishnevsky D., Lisitsa V. Lebedev scheme for numerical simulation of wave propagation in 3D anisotropic elasticity. Geophysical Prospecting 58, 2010, p. 619-635. DOI 10.1111/j.1365-2478.2009.00862.x
12. William Gropp, Steven Huss-Lederman et al. MPI - The complete reference. Cambridge, Massachusetts, MIT Press, 1998, vol. 2, 205 p.
13. Glinskii B. M., Karavaev D. A., Kovalevsky V. V., Martynov V. N. Numerical modeling and experimental research of the «Karabetov Mountain» mud volcano by vibroseismic methods. Num. Meth. Prog., 2010, vol. 11 (1), p. 95-104.
14. Collino F., Tsogka C. Application of the Perfectly Matched Absorbing Layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media. Geophysics 66, 2001, p. 294-307. DOI 10.1190/1.1444908
References
1. Liseikin V. D. Grid generation method. 2nd ed. Berlin, Springer, 2010. ISBN 978-90-4812912-6
2. Vaseva I. A., Liseikin V. D., Likhanova Yu. V., Morokov Yu. N. An elliptic method for construction of adaptive spatial grids. Russ. J. Numer. Anal. Math. Modelling, 2009, vol. 1, p. 65-78. DOI 10.1515/RJNAMM.2009.006
3. Liseikin V. Differential meshes. Theory and applications. Novosibirsk, SB RAS Publ., 2014, 254 p. (in Russ.)
4. Appelo D., Petersson N. A. A stable finite difference method for the elastic wave equation on complex geometries with free surfaces. Communications in computational physics 5, 2007, no. 1, p. 84-107.
5. Jose M. Carcione. The wave equation in generalized coordinates. Geophysics 59, 1994, no. 12, p. 1911-1919. DOI 10.1190/1.1443578
6. Titov P. An algorithm and a program for simulation of 2D-wave fields in areas with a curved free surface. In: Proceedings of "Scientific service on the Internet - 2014" conference. Novorossiysk, Abrau-Dyurso, 2014, p. 446-455. (in Russ.)
7. Titov P. Modeling of elastic waves in media with a complex free surface topography. Vestnik of NSU. Series: Information Technologies, 2018, vol. 16, no. 4, p. 153-166. (in Russ.)
8. Titov P. A technology of modeling of elastic waves propagation in media with complex 3D geometry of the surface with the aid of high performance cluster. Russian Supercomputing Days, Proc, 2016, p. 1020-1031.
9. Titov P. The Simulation of 3D Wave Fields in Complex Topography Media. In: Russian Supercomputing Days, CCIS 1129, 2019, p. 451-462. DOI 10.1007/978-3-030-36592-9_37
10. Vishnevsky D., Lisitsa V., Tcheverda V., Reshetova G. Numerical study of the interface errors of finite-difference simulations of seismic waves. Geophysics, 2014, vo. 79 (4), p. 219232. DOI 10.1190/geo2013-0299.1
11. Vishnevsky D., Lisitsa V. Lebedev scheme for numerical simulation of wave propagation in 3D anisotropic elasticity. Geophysical Prospecting 58, 2010, p. 619-635. DOI 10.1111/j.1365-2478.2009.00862.x
12. William Gropp, Steven Huss-Lederman et al. MPI - The complete reference. Cambridge, Massachusetts, MIT Press, 1998, vol. 2, 205 p.
13. Glinskii B. M., Karavaev D. A., Kovalevsky V. V., Martynov V. N. Numerical modeling and experimental research of the «Karabetov Mountain» mud volcano by vibroseismic methods. Num. Meth. Prog., 2010, vol. 11 (1), p. 95-104.
14. Collino F., Tsogka C. Application of the Perfectly Matched Absorbing Layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media. Geophysics 66, 2001, p. 294-307. DOI 10.1190/1.1444908
Материал поступил в редколлегию Received 09.08.2020
Сведения об авторе
Титов Павел Андреевич, младший научный сотрудник лаборатории Суперкомпьютерного моделирования, Институт вычислительной математики и математической геофизики СО РАН (Новосибирск, Россия)
Information about the Author
Pavel A. Titov, Junior Researcher, Supercomputing Modeling Laboratory, Institute of Computational mathematics and Mathematical Geophysics SB RAS (Novosibirsk, Russian Federation)
