автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численное моделирование сейсмических процессов на высокопроизводительных вычислительных системах

005004193

На правах рукописи

ХОХЛОВ Николай Игоревич

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ СЕЙСМИЧЕСКИХ ПРОЦЕССОВ НА ВЫСОКОПРОИЗВОДИТЕЛЬНЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ

Специальность 05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

- 1 ДЕК 2011

МОСКВА-2011

005004193

Работа выполнена на кафедре информатики Московского физико-технического института (государственного университета)

Научный руководитель:

доктор физико-математических наук, профессор Петров Игорь Борисович

Официальные оппоненты:

доктор физико-математических наук, профессор Трошкин Олег Валентинович

кандидат технических наук, ст. н. с. Нигметов Геннадий Максимович

Ведущая организация:

Институт прикладной математики им. М. В. Келдыша РАН

Защита состоится «/Л » 2011 г. в ^ часов на

заседании диссертационного совета Д 212.156.05 при Московском физико-техническом институте (государственном университете) по адресу: 141700, Московская обл., г. Долгопрудный, Институтский пер., д. 9, ауд. 903 КПМ.

С диссертацией можно ознакомиться в библиотеке МФТИ. Автореферат разослан «/%> 2011 г.

Ученый секретарь диссертационного совета

Федько О. С.

Общая характеристика работы

Актуальность темы

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

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

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

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

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

Цели работы

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

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

Научная новизна

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

2. Проведен сравнительный анализ известных лимитеров для ТУБ схем второго порядка точности и схем до 4-го порядка точности. На основе полученных результатов выбрана схема, позволяющая получать достаточно точные численные решения в областях больших градиентов.

3. Получено и исследованно поведение поверхностных волн Лява и Рэлея при численном решении динамической пространственной задачи теории упругости.

4. Проведено моделирование процесса воздействия фронта упругой волны, инициированной землетрясением, через наземное сооружение, находящееся во вмещающем массиве. Определены области разрушения здания.

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

6. Предложен способ распараллеливания алгоритма для работы в распределенной кластерной среде на основе технологии MPI. Для распараллеливания работы программы в системах с общей памятью использована технология ОрепМР.

7. Проведена оптимизация алгоритма решения систем гиперболических уравнений на регулярных сетках для работы на центральных процессорах с эффективным использованием кеша и потоковых инструкций S SE.

8. Разработан механизм задания геометрии и входных параметров среды, использующий свободный пакет для пре- и постпроцессинга Gmsh.

Практическая ценность

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

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

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

Работа выполнена при поддержке ряда государственных и коммерческих грантов:

• Программа (мероприятие): федеральная целевая программа «Научные и научно-педагогические кадры инновационной России» на

2009-2013 гг., в рамках реализации мероприятия № 1.2.1 «Проведение научных исследований научными группами под руководством докторов наук». Проект: «Разработка вычислительных технологий для моделирования пространственных динамических процессов в проблеме сейсморазведки на высокопроизводительных ЭВМ»;

• Федеральное государственное унитарное предприятие «Российский Федеральный Ядерный Центр - Всероссийский научно-исследовательский институт экспериментальной физики (ФГУП «РФЯЦ-ВНИИЭФ»)». НИР5. «Разработка физико-математических моделей, алгоритмов и эффективных методов решения задач механики сплошных сред на супер-ЭВМ»;

• Грант РФФИ 11-01-12011-офи-м-2011. Разработка численных методов для решения задач геомеханики и сейсморазведки на многопроцессорных вычислительных системах, 2011-2012 гг.;

• Грант РФФИ 10-01-92654-ИНД_а. Математическое моделирование сложных задач на высокопроизводительных вычислительных системах.

2010-2011 гг;

• Грант РФФИ 11-08-01240-а. Моделирование ударных, взрывных и сейсмических воздействий на сооружения атомной энергетики. 2011 -2012 гг;

• Small or medium-scale focused research project (STREP) proposal ICT EU-Russia Coordinated Call. FP7-2011-EU-Russia. 2011-2012 y.

Публикации

Научные результаты диссертации опубликованы в 23 работах, из которых две [8, 15] - в изданиях, рекомендованных ВАК для публикации основных результатов диссертации.

Апробация

Результаты работы были доложены, обсуждены и получили одобрение специалистов на следующих научных конференциях:

• Научные конференции Московского физико-технического института - Всероссийские молодёжные научные конференции с международным участием «Проблемы фундаментальных и прикладных, естественных и технических наук в современном информационном обществе» (МФТИ, Долгопрудный, 2006 - 2009);

• XVI международная научно-практическая конференция «Комплексная безопасность 2011 г.» (МЧС, Москва, 2011);

• VI международная конференция по неравновесным процессам в соплах и струях (Санкт-Петербург, 2006);

• 20-th International Conference On Transport Theory (ICTT-20) (Обнинск, 2007);

• 26-th International Symphosium on Rarefied Gas Dynamics (Kyoto University, Japan, 2008);

• XII международный семинар «Супервычисления и математическое моделирование» (РФЯЦ - ВНИИЭФ, Саров, 2011);

• APOS-EU/APOS-RU Project Workshop "Modelling of seismic problems" (ИПМ РАН, Москва, 2011);

• NPP safety and personal training. XII International Conference. (Обнинск, 2011);

• Российско-индийский семинар «Новые достижения математического моделирования» (ИАП РАН, Москва, 2011).

Результаты работы были доложены, обсуждены и получили одобрение специалистов на научных семинарах в следующих организациях:

• Всероссийский научно-исследовательский институт по проблемам гражданской обороны и чрезвычайных ситуаций МЧС России (федеральный центр) (Москва, 2011);

• Институт прикладной математики им. М. В. Келдыша РАН (Москва, 2011);

• Российский федеральный ядерный центр - Всероссийский научно-исследовательский институт экспериментальной физики (Саров, 2011).;

• Московский государственный университет путей сообщения (Москва, 2011).

Структура и объем диссертации

Диссертация состоит из введения, пяти глав и заключения. Общий объем диссертации составляет 110 страниц. Список литературы содержит 70 публикаций.

Содержание работы

Введение

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

Глава 1

В первой главе описаны основные цели диссертации, приведены постановки задач, решенных в работе.

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

Вторая задача заключалась в исследовании поверхностных волн Лява и Рэлея, изучении их характеристик и сравнении полученных результатов с теоретическими. Задача представляется актуальной вследствие большого значения данных волн в различных сейсмических

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

Третья задача представляла собой моделирование распространения волн, инициированных землетрясением, от гипоцентра до земной поверхности. Для моделирования очага землетрясения была выбрана сдвиговая модель. Такая модель физически соответствует ситуации, когда в земной коре существует разлом, по которому и происходит подвижка при землетрясении. Выбор величины модуля скорости был произведён путём сопоставления данных моделирования и реальных экспериментальных данных. Было проведено моделирование землетрясения во вмещающем массиве слоистой среды, плотность каждого слоя принималась равной 2500 кг/м3. В расчете задавались 4 слоя, скорость звука в слоях увеличивалась с глубиной слоя. В этой задаче основной интерес представляют волновые процессы в массиве земной породы, которые инициируются землетрясением, а также модель для задания гипоцентра землетрясения и прохождение волн через слоистые структуры.

Глава 2

Во второй главе описаны определяющие уравнения и приведена математическая модель деформируемого твердого тела.

Определяющие уравнения

Рассмотрим основные уравнения линейной динамической теории упругости. Это уравнения движения и реологические соотношения:

"В=7Т- (1) Т = А(У-й)1 + Му®^ + и®у)'

здесь р - плотность среды, й - вектор скорости смещения, Т - тензор напряжений Коши, V - градиент по пространственным координатам, Л и ц - упругие постоянные Ляме, I - единичный тензор, ® - оператор тензорного произведения. Второе уравнение представляет собой закон Гуна, продифференцированный по времени.

Составим вектор искомых функций, состоящий из 9-ти

компонент:

и = ,

здесь о1 - компоненты вектора скорости смещения, ст. - компоненты

симметричного тензора напряжений Т. В этом случае допускается запись системы уравнений (1) динамики деформируемого твердого тела в матричном виде:

Й1 Эч

— = 2А—<- (2)

а« ОХ.

где Ау - матрицы размера 9x9.

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

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

Системы гиперболических уравнений

Для численного моделирования задач динамики деформируемого твердого тела широко применяется класс гибридных сеточно-характеристический метод, что обеспечивает моность численных решений. После замены переменных п = \Уи каждая из

систем (2) распадается на девять независимых скалярных уравнений переноса (индекс j далее опускается там, где это возможно):

Зп дп

— + L— = 0.

dt дх

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

u = W п .

Разностные схемы

Далее в работе описаны реализованные конечно-разностные численные методы, сформулированные для гиперболических систем уравнений с одной пространственной переменной: схемы Куранта-Изаксона-Риса, Лакса-Вендорфа, Бима-Уорминга, Фромма, симметричная схема 4-го порядка точности, гибридная схема Федоренко 2-го порядка точности и TVD-схема 2-го порядка точности с большим набором ограничителей. В работе приведен сравнительный анализ следующих ограничителей: superbee, MC, Minmod, Koren, Osher, Smart, Sweby, UMIST, van Albada 2, van Leer.

Среди реализованных схем по результатам проведенных тестов наилучшее качество решения обеспечивает TVD-схема с ограничителем superbee:

здесь /т+1/2 и /т_,., - антидиффузионые потоки:

/.♦w: = «'. +0.5ф(гяХ1-егХи'^ - О,

-о

Функция ф(г ) называется ограничителем или лимитером. Параметр гт вычисляется по формуле:

и -и ,

г = —1--

»«i

Основной ограничитель, использованный во всех расчетах, - superbee, он имеет следующий вид:

ф{г) = max[0,min(2r,l), min(r,2)].

На рис. 1 показан сравнительный график точного решения и решения, полученного при использовании различных разностных схем.

I ТНео'е'|са1 ■■

С1А -| ТУО зирегЬее , НуЬгй (гесогепКо) 4^1 огйег -

II

ЬЬ ОУ ^

Рис 1: Сравнение различных разностных схем

На данном графике представлены результаты для схемы первого порядка точности Куранта-Изаксона-Риса, ТУБ-схемы второго порядка точности с ограничителем зирегЬее, гибридной схемы второго порядка точности Федоренко и схемы четвертого порядка точности. Как видно наилучшее сходство с решением дает ТУБ-схема с ограничителем зирегЬее. В работе представлено сравнение 19 вариации разностных схем и лимитеров, и, в большинстве случаев, ограничитель вирегЬее дал лучшее решение на тестовых расчетах.

Криволинейные сетки

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

Расщепление по координатам производиться в этом случае не в физическом, а в преобразованном пространстве. Исходная система уравнений (2) преобразуется к виду:

Теперь для каждого шага расщепления для направления матрица L будет иметь вид:

L = diag (с/., , cj., -с2/., с2/,, -с2Л, 0,0,0}, здесь I =|ny |= ^(v;)2+(vJ2f+(v'3)\ a Матрицы

преобразования W и W"1 также изменяются.

Для определения коэффициентов nJ вычисляется матрица, обратная к матрице Якоби, ее необходимо вычислять в каждой ячейке. Расчетная область

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

системы координат: " = Любой из блоков Qt характеризуется

своим набором параметров среды, граничных и контактных условий. В каждом блоке Q, строятся прямоугольные сетки с постоянным внутри

блока шагом h4 = {hkl,hkl,hj, при этом в соседних блоках, при

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

В четвертой главе описан комплекс программ, особенности его реализации, оптимизация алгоритма для работы на центральных процессорах, технология распараллеливания с использованием библиотек MPI и ОрепМР. Также рассмотрены способы хранения данных в памяти в параллельной версии программы с целью оптимизации времени расчета.

Оптимизация работы с кешем процессора и использование

SSE

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

Несмотря на то, что в расчетах используется явная разностная схема, нет необходимости хранить расчетную сетку для двух временных слоев. За счет использования метода расщепления по пространству достаточно хранить всего лишь несколько расчетных улов сетки для перехода на следующий временной шаг. Однако из-за особенности хранения сетки в языках C/C++ (в виде одномерного массива) при такой схеме возникает много промахов кеша при расчетах. В работе описан алгоритм оптимизации работы с кешем за счет несколько большего потребления оперативной памяти, однако это позволило уменьшить

время работы тестовой программы до 5 раз.

220

200

180

160

§ 140 с

& 120

I 100

S 80 60 40 20 0

0 1000 2000 3000 4000 5000 6000 7000

Grid size

Рис 2: График зависимости числа тактов на перенос одного узла от размера сетки

Далее описана технология использования SSE-инструкций в применяемом алгоритме для ускорения вычислений. За счет некоторого изменения структуры хранения данных на основе чисел типа float удалось эффективно использовать потоковые инструкции процессоров, что дало дополнительное ускорение до двух раз. На рис. 2 показаны графики зависимости числа тактов процессора, которые процессор тратит на пересчет одного узла, в зависимости от размера сетки по одной из осей. В тестовых расчетах использовалось двумерное

Core2 Quad Q8200 (Not optimized) - Core2 Quad Q8200 (Optimized)--- Cor©2 Quad Q8200 (Optimized +■ SSE)

... ^ f. -

'.^Vi"''

. , ,-------¡-""""j ; I 1 f

■X..............-i.............. !! ^

уравнение акустики. Сходная картина получается для различных версий компиляторов (gcc, icc) и процессоров.

В расчетной программе для уравнений теории упругости была использована только часть описанной выше оптимизации, и это дало ускорение в 3 раза на схеме первого порядка точности и в 2 раза на схеме второго порядка точности.

Распараллеливание

В программном комплексе реализованы две технологии распараллеливания. Для расчетов, ведущихся в общей памяти, используется технология ОрепМР, на распределенных кластерных системах - технология MPI. За счет использования технологии MPI удалось значительно увеличить скорость расчета задач и увеличить размеры используемых сеток. Алгоритм построен таким образом, что не имеет ограничений по максимальному числу расчетных узлов. Все тестовые расчеты проводились на сетках до миллиарда узлов. На рис. 3 показан график ускорения расчета при росте числа расчетных ядер на основе технологии MPI.

50 45

35 30

а.

! 25 а. <п

20 15 10 5 О

Number proceses

Рис 3: График зависимости ускорения расчета в зависимости от числа ядер

В тесте использовалась сетка из 64 млн узлов, расчеты велись на кластере МФТИ-60.

Комплекс программ

Далее в работе описана структура реализованного комплекса программ. В качестве инструментов для пре- и постпроцессинга используются свободные программные пакеты Gmsh, ParaView и Mayavi

о 10 20 30 40 50 60 70

с доработками в виде набора скриптов и подпрограмм. Данные сохраняются в широко используемый формат визуализации научных данных УТК. Сам программный комплекс имеет широкий набор параметров входных данных, расчетных геометрий, граничных условий и сохранения результатов. Более подробно возможности и особенности программного комплекса описаны в диссертации.

Глава 5

В пятой главе приведены результаты численного моделирования для нескольких задач, анализ полученных результатов и сравнение ряда результатов с аналитическими.

Рис 4: Результаты расчета для поперечной волны: а) I = 0,0013 с. Ь) 1 - 0,0039 с, с) I —

0,0079 с

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

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

Также моделировалось распространение волн Лява и Рэлея. Для волн Лява расчетная область была разделена на два слоя с различными параметрами среды и между ними была выделена контактная граница. На рис. 5 показан результат моделирования распространения волны Лява в различные моменты времени. Градациями серого цвета отображена горизонтальная составляющая скорости волны, перпендикулярная направлению распространения волны. Хорошо видно, что волна имеет горизонтальную поляризацию, очень быстро затухает в нижнем полупространстве, по мере распространения волны ее форма меняется.

Рис 5: Распространение волны Лява Распространение волны Рэлея представлено на рис. 6. На графике показано распределение модуля скорости среды. Расчетная скорость волны показала совпадение с теоретической.

оретиче'

Рис 6: Распространение волны Рэлея Следующий набор расчетов связан с моделированием распространение волн от гипоцентра землетрясения до дневной поверхности земли. В работе используется модель задания возмущений в гипоцентре землетрясений: возмущения задается путем сдвига

соседних областей в противоположных направлениях. Область гипоцентра рассматривалась в виде параллелепипа во всех расчетах. Это дает хорошее качественное сходство с реальными данными.

ЕаКИчиаке 1

ЕаППЧиаке 2

Ва»

(¡1

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

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

Рис 8: Распространение возмущений от точечного взрыва и землетрясений с двумя различными параметрами Моделирование проводилось в многослойной модели земной коры, на сейсмограммах хорошо заметны отражения от различных слоев. На рис. 8 показаны распределения скоростей для точечного взрыва и двух типов землетрясений в расчетной области в один и тот же момент времени.

Также в работе представлено еще несколько тестовых расчетов.

В заключении приводятся основные результаты и выводы диссертации.

Основные результаты и выводы диссертации

1. Создан комплекс программ для решения пространственных задач динамики деформируемого твердого тела. Создана программная среда для пре- и постпроцессинга. Созданный комплекс позволяет моделировать широкий круг задач сейсмики и сейсмостойкости.

2. Численно получено и исследовано поведение поверхностных волн Лява и Рэлея. Показана соответствие полученных расчетных результатов с теоретическими.

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

4. Исследованы особенности распространения волн, инициированных землетрясением, от гипоцентра до земной поверхности во многослойных геологических средах.

5. Разработаны и реализованы гибридные сеточно-характеристические схемы, с различными ограничителями, от первого до четвертого порядка точности для численного решения пространственных задач сейсмостойкости.

6. Выполнено распараллеливание комплекса расчетных программ на основании технологий MPI и OpenMP. ~

7. Проведена оптимизация расчетного алгоритма для работы на современных центральных процессорах, что позволило повысить скорость работы алгоритма.

Список работ по теме диссертации

1. Golubev V. /., Kvasov I. Е., Petrov I. В., Khokhlov N. I. Numerical modeling of dynamic processes occurring under the influence of earthquakes on the power structures in high-performace computers // NPP safety and personal training. XII International Conference Anstracts - Obninsk, 2011. P. 109-110.

2. Golubev V. I., Petrov I. В., Khokhlov N. I. Simulation of seismic ground structures using the grid-characteristic method // Russian-indian workshop. Book of abstracts. Advansced Computational Modeling and Simulations -Moscow, 2011. -P. 7.

3. Khokhlov N. I., Kloss Yu. Yu., Shurigin B. A., Tcheremissine F. G. Application of MPI-technology for solving the Boltzmann equation // 20-th International Conference on Transport Theory - Obninsk, 2007. - P. 189.

4. Khokhlov N. I., Kloss Yu. Yu., Shurigin B. A., Tcheremissine F. G. Simulation of the temperature driven micro pump by solving the Boltzmann equation // 26-th International Symphosium on Rarefied Gas Dynamics -2008.-P. 115.

5. Голубев В. И., Квасов И. Е, Хохлов Н. И., Петров И. Б. Численное моделирование волн Рэлея и Лява сеточно-характеристическим методом. // Сборник материалов XVI международной научно-практической конференции "Комплексная безопасность 2011 г." - 2011.

- С. 25-27.

6. Голубев В. И., Квасов И. Е., Хохлов Н. И., Петров И. Б. Численное моделирование сейсмостойкости жилых и промышленных наземных сооружений. // Сборник материалов XVI международной научно-практической конференции "Комплексная безопасность 2011 г.", - 2011.

- С. 27-30.

7. Дербакова Е. П., Кчосс Ю. Ю„ Хохлов Н. И., Федотов В. Ю„ Шурыгин Б. А. Паралелльные алгоритмы численных схем решения уравнения Больцмана на основе технологии MPI // Электронный журнал "Исследовано в России" - 2007. - № 54 - С. 581-588.

8. Клосс Ю. Ю., Черемисин Ф. Г., Хохлов Н. К, Шурыгин Б. А. Программно-моделирующая среда для исследования течений газа в микро- и наноструктурах на основе решения уравнения Больцмана // Атомная энергия - 2008. - № 105 Т. 4 - С. 211-213.

9. Клосс Ю. Ю., Дербакова Е. П., Хохлов Н. И., Федотов В. Ю. Тестирование сеточных схем с параллельной реализацией для уравнения Больцмана на кластерных системах. // Современные проблемы фундаментальных и прикладных наук: Труды 50-й научной конференции. / Моск. физ. - техн. ин-т. - М. - Долгопрудный, 2007 - С. 9192.

10.Клосс Ю. 10., Хохлов Н. И. Технологии распределенных вычислений для сеточных схем численного решения уравнения Больцмана на основе библиотеки MPI. // Современные проблемы фундаментальных и прикладных наук: Труды 50-й научной конференции. / Моск. физ. - техн. ин-т. - М. - Долгопрудный, 2007 - С. 8991.

11. Клосс Ю. Ю., Черемисин Ф. Г., Хохлов Н. И., Гришина В. Г., Сакмаров А. В. Реализация консервативной схемы первого порядка в распределенной вычислительной среде на основе стандарта MPI // Сборник 6-й Международной конференции по неравновесным процессам в соплах и струях - Санкт-Петербург, 2006 - С. 146-149.

12. Класс Ю. Ю„ Черемисин Ф. Г., Хохлов Н. И., Шурыгин Б. А. Компьютерное моделирование вакуумного микронасоса без движущихся частей на основе уравнения Больцмана. // Современные проблемы фундаментальных и прикладных наук: Труды 50-й научной конференции. / Моск. физ. - техн. ин-т. - М. - Долгопрудный, 2007 - С. 9593. „ ,

13.Клосс Ю. Ю., Черемисин Ф. Г., Хохлов Н. И., Шурыгин Б. А. Разработка численных схем решения кинетического уравнения в кластерных средах на основе технологии MPI. // Информационные процессы, - 2007. - №7 Т. 4 - С. 425-431.

14. Клосс Ю. Ю., Черемисин Ф. Г., Хохлов Н. К, Шурыгин Б. А., Шувалов П. В. Моделирование и анализ газокинетических процессов в микро- и наноструктурах. Книга из серии из 12 книг по направлению "Физико-технические проблемы ядерной энергетики" под редакцией академика РАН H. Н. Пономарева-Степного. - М.:РНЦ Курчатовский

институт, 2008. - 48 с.

15. Петров И. Б., Хохлов Н. И. Моделирование сейсмических явлений сеточно-характеристическим методом. // Труды МФТИ, - 2011. - № 3 Т. 3-С. 159-167.

16. Петров И. Б., Хохлов Н. И. Сравнение TVD лимитеров для численного решения уравнений динамики деформируемого твердого тела сеточно-характеристическим методом. // Математические модели и задачи управления.: Сборник научных трудов/ Моск. физ.-тех. ин-т., -

M., 2011.-С. 104-111.

17. Субботина А. Ю., Хохлов Н. И. Реализация клеточных автоматов "игра "Жизнь"" и клеточного автомата Кохомото-Ооно с применением технологии MPI. // Компьютерные исследования и моделирование, -2010.- № 2 Т. 3 - С. 319-322.

18. Хохлов Н. И. Распараллеливание метода решения задач моделирования волновых процессов в деформируемом твердом теле используя технологию MPI. // Математические модели и задачи управления.: Сборник научных трудов/ Моск. физ.-тех. ин-т., - М., 2011. -С. 112-116.

19.Хохлов Н. И. Тестирование латентности и скорости обмена данными библиотеки МР1 при использовании сети Муппй. // Модели и методы обработки информации Сборник научных трудов/ Моск. физ.-тех. ин-т., -М., 2009. - С. 107-115.

20.Хохлов Н. И., Гришина В. Г., Федотов В. Ю. Применение численных сеточных схем решения уравнения Больцмана для моделирования и анализа течений газа в микроустройствах и наноструктурах. // Современные проблемы фундаментальных и прикладных наук: Труды 51-й научной конференции. / Моск. физ. - техн. ин-т. - М. - Долгопрудный, 2008 - С. 62-63.

21. Хохлов Н. И., Дербакова Е. П., Клосс Ю. Ю. Компьютерные модели на основе уравнения Больцмана для параметрического исследования микронасоса без движущихся частей. // Современные проблемы фундаментальных и прикладных наук: Труды 51-й научной конференции. / Моск. физ. - техн. ин-т. - М. - Долгопрудный, 2008 - С. 6466.

22. Черемисин Ф. Г., Шурыгин Б. А., Хохлов Н. И. Моделирование и анализ течений разреженного газа сквозь периодическую наноструктуру. // Современные проблемы фундаментальных и прикладных наук: Труды 49-й научной конференции. / Моск. физ. - техн. ин-т. - М. - Долгопрудный, 2006 - С. 89-90.

23.Шурыгин Б. А., Хохлов Н. И. Создание интерактивной оболочки для визуализации и анализа газокинетических процессов тепломассопереноса в микро и наноструктурах. // Современные проблемы фундаментальных и прикладных наук: Труды 49-й научной конференции. / Моск. физ. - техн. ин-т. - М. - Долгопрудный, 2006 - С. 9192.

В работах с соавторами лично соискателем выполнено следующее:

1. Создан комплекс программ для численного решения пространственных динамических задач механики деформируемого твердого тела с учетом неоднородностей среды, наличием подземных и наземных сооружений.

2. Проведено моделирование распространения поверхностных сейсмических волн Лява и Рэлея, сравнение полученных результатов с теорией.

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

4. Проведено моделирование динамических процессов в многослойной коре земли, инициированных землетрясением, от гипоцентра землетрясения до земной поверхности.

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

6. Проведен сравнительный анализ лимитеров для TVD-схем (второго порядка точности) и схем до четвертого порядка точности применительно к данному типу задач.

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

8. Реализовано распараллеливание комплекса программ для работы на многоядерных многопроцессорных машинах, в распределенной кластерной среде.

9. Алгоритм оптимизирован для эффективной работы на больших сетках (до миллиарда узлов), проведена адаптация расчетного модуля для эффективного использования кеша процессора и ускорения расчетов с использованием потоковых инструкций SSE.

Подписано в печать: 07.11.2011 Объем: 1,5 усл.п.л. Тираж: 100 экз. Заказ № 539 Отпечатано в типографии «Реглет» 119526, г. Москва, Страстной бульвар, д. 6,стр. (495) 978-43-34; \v\vw.reglct.ru

Введение

Общая характеристика работы.

Методы моделирования сейсмических процессов.

Механизмы очага землетрясения.

Численные методы решения уравнений гиперболического типа

Глава 1. Цели диссертации и решаемые задачи

Цели диссертационной работы.

Прохождение сейсмических волн через наземные сооружения

Сейсмические волны

Распространение поверхностных волн Лява и Рэлея.

Распространение волн от гипоцентра землетрясения

Глава 2. Определяющие уравнения

Определяющие уравнения

Выбор системы координат

Обобщенная запись

Случай параллелепипедной сетки

Глава 3. Численные методы, используемые в работе.

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

Системы линейных уравнений гиперболического типа

Вид матриц Л, Г2, ГІ

Метод расщепления по пространственным координатам

Сеточно-характеристический метод

Криволинейные расчетные сетки

Разностные схемы, исследуемые в работе.

Разностные схемы для одномерного уравнения переноса

TVD-разностные схемы второго порядка точности

Результаты тестирования разностных схем.

Глава 4. Программный комплекс

Расчетная область.

Особенности хранения расчетных сеток.

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

Описание.

Оптимизация работы с памятью.

Оптимизация sse.

Выводы.

Распараллеливание

Распараллеливание на основе технологии MPI.

Алгоритм деления расчетной области.

Обмен граничными значениями

Результаты тестирования.

Распараллеливание на основе технологии ОрепМР

Выводы.

Программный комплекс.

Архитектура.

Вычислительный модуль

Препроцессинг.

Постпроцессинг

Глава 5. Результаты численного моделирования ряда задач

Волны Лява

Описание.

Результаты моделирования.

Волны Рэлея.

Описание.

Результаты моделирования.

Прохождение волн через наземные объекты

Описание.'.

Результаты расчетов

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

Описание.

Параметры расчетов

Результаты расчета

Прохождение волн через водонасыщенные грунты

Описание.

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

Интерполяция расчетных данных

Изменение реологических свойств водонасыщенных грунтов. Уравнения Гассмана

Результаты моделирования.

Общая характеристика работы

Актуальность темы

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

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

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

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

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

Научная новизна

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

• Проведен сравнительный анализ известных лимитеров для ТУБ схем второго порядка точности и схем до 4-го порядка точности. На основе полученных результатов выбрана схема, позволяющая получать достаточно точные численные решения в областях больших градиентов.

• Получено и исследованно поведение поверхностных волн Лява и Рэлея при численном решении динамической пространственной задачи теории упругости.

• Проведено моделирование процесса воздействия фронта упругой волны, инициированной землетрясением, через наземное сооружение, находящееся во вмещающем массиве. Определены области разрушения здания.

• Найдено численное решение задачи распространения волн, инициированных динамическим возмущением, в земной коре через водонасыщен-ные грунты.

• Предложен способ распараллеливания алгоритма для работы в распределенной кластерной среде на основе технологии MPI. Для распараллеливания работы программы в системах с общей памятью использована технология ОрепМР.

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

• Разработан механизм задания геометрии и входных параметров среды, использующий свободный пакет для пре- и постпроцессинга Gmsh.

Практическая ценность

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

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

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

Работа выполнена при поддержке ряда государственных и коммерческих грантов:

• Программа (мероприятие): федеральная целевая программа «Научные и научно-педагогические кадры инновационной России» на 2009-2013 гг., в рамках реализации мероприятия № 1.2.1 «Проведение научных исследований научными группами под руководством докторов наук». Проект: «Разработка вычислительных технологий для моделирования пространственных динамических процессов в проблеме сейсморазведки на высокопроизводительных ЭВМ»;

• Федеральное государственное унитарное предприятие «Российский Федеральный Ядерный Центр - Всероссийский научно-исследовательский институт экспериментальной физики (ФГУП «РФЯЦ-ВНИИЭФ»)». НИР5 «Разработка физико-математических моделей, алгоритмов и эффективных методов решения задач механики сплошных сред на супер-ЭВМ»;

• Грант РФФИ 11-01-12011-офи-м-2011. Разработка численных методов для решения задач геомеханики и сейсморазведки на многопроцессорных вычислительных системах, 2011-2012 гг.;

• Грант РФФИ 10-01-92654-ИНДа. Математическое моделирование сложных задач на высокопроизводительных вычислительных системах. 2010

- 2011 гг;

• Грант РФФИ 11-08-01240-а. Моделирование ударных, взрывных и сейсмических воздействий на сооружения атомной энергетики. 2011 - 2012 гг;

• Small or medium-scale focused research project (STREP) proposal ICT EU-Russia Coordinated Call. FP7-2011-EU-Russia. 2011-2012 y.

Публикации

Научные результаты диссертации опубликованы в 23 работах, из которых две [48, 51] - в изданиях, рекомендованных ВАК для публикации основных результатов диссертации.

В работах с соавторами лично соискателем выполнено следующее:

• Создан комплекс программ для численного решения пространственных динамических задач механики деформируемого твердого тела с учетом иеоднородностей среды, наличием подземных и наземных сооружений.

• Проведено моделирование распространения поверхностных сейсмических волн Лява и Рэлея, сравнение полученных результатов с теорией.

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

• Проведено моделирование динамических процессов в многослойной коре земли, инициированных землетрясением, от гипоцентра землетрясения до земной поверхности.

• Получено численное решение задачи распространения волн, инициированных динамическим возмущением в земной коре, через водонасыщен-ные грунты.

• Проведен сравнительный анализ лимитеров для TVD-схем (второго порядка точности) и схем до четвертого порядка точности применительно к данному типу задач.

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

• Реализовано распараллеливание комплекса программ для работы на многоядерных многопроцессорных машинах, в распределенной кластерной среде.

• Алгоритм оптимизирован для эффективной работы на больших сетках (до миллиарда узлов), проведена адаптация расчетного модуля для эффективного использования кеша процессора и ускорения расчетов с использованием потоковых инструкций SSE.

Апробация

Результаты работы были доложены, обсуждены и получили одобрение специалистов на следующих научных конференциях:

Научные конференции Московского физико-технического института -Всероссийские молодёжные научные конференции с международным участием «Проблемы фундаментальных и прикладных, естественных и технических наук в современном информационном обществе» (МФТИ, Долгопрудный, 2006 - 2009);

XVI международная научно-практическая конференция «Комплексная безопасность 2011 г.» (МЧС, Москва, 2011);

VI международная конференция по неравновесным процессам в соплах и струях (Санкт-Петербург, 2006);

20-th International Conference On Transport Theory (ICTT-20) (Обнинск, 2007);

26-th International Symphosium on Rarefied Gas Dynamics (Kyoto University, Japan, 2008);

XII международный семинар «Супервычисления и математическое моделирование» (РФЯЦ - ВНИИЭФ, Саров, 2011);

APOS-EU/APOS-RU Project Workshop "Modelling of seismic problems" (ИШ РАН, Москва, 2011);

NPP safety and personal training. XII International Conference. (Обнинск, 2011);

Российско-индийский семинар «Новые достижения математического моделирования» (ИАП РАН, Москва, 2011).

Результаты работы были доложены, обсуждены и получили одобрение специалистов на научных семинарах в следующих организациях:

• Всероссийский научно-исследовательский институт по проблемам гражданской обороны и чрезвычайных ситуаций МЧС России (федеральный центр) (Москва, 2011);

• Институт прикладной математики им. М. В. Келдыша РАН (Москва, 2011);

• Российский федеральный ядерный центр - Всероссийский научно-исследовательский институт экспериментальной физики (Саров, 2011).;

• Московский государственный университет путей сообщения (Москва, 2011).

Структура и объем диссертации

Диссертация состоит из введения, пяти глав и заключения. Общий объем диссертации составляет 110 страниц. Список литературы содержит 70 публикаций.

Выводы

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

При распараллеливании удалось добиться эффективности до 76% при увеличении числа ядер в 64 раза. Применять технологию ОрепМР при расчетах на данном программном комплексе стоит в тех случаях, когда эффективность MPI уже достигла своего минимума. Тогда дальнейшее ускорение может быть достигнуто за счет технологии ОрепМР, при использовании данной технологии внутри одного расчетного узла на гибридных вычислительных машинах.

Применение технологий MPI и ОрепМР позволило существенно уменьшить время расчета и увеличить размеры используемых расчетных сеток.

Программный комплекс Архитектура

Одним из результатов данной работы является программный комплекс, который позволяет численно решать задачи распространения волн, инициируемых различными динамическими нагрузками на трехмерные объекты со сложной внутренней структурой. Разработанная параллельная версия комплекса программ позволяет существенно сократить время вычислений, а также увеличить разрешение моделей. Программный комплекс применим также для решения задач ЗБ-сейсмики и сейсморазведки.

На рис. 4.9 изображена архитектура программного комплекса, он состоит из:

• вычислительного модуля, основанного на сеточно-характеристическом

• расчет на блочной геометрии из комбинации сеток, в местах контактов сеток они должны быть согласованы;

• включает в себя монотонные сеточно-характеристические схемы 1-4 порядка точности;

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

Архитектура ядра вычислительного комплекса разрабатывалась с учетом дальнейших расширений возможностей программного комплекса. При написании программ существует эмпирическое правило 80-20, которое говорит о том, что 80% исходного кода выполняется 20% времени и, соответственно, наоборот. Процесс разработки и отладки был значительно ускорен и упрощен, за счет написания частей кода, которые не критичны к скорости выполнения или выполняются один раз за расчет, на скриптовом языке высокого уровня Python.

Вначале все важные модули, критичные к скорости выполнения, были написаны и отлажены на языке С/С-Н-, к ним относится расчетная сетка, граничные и контактные условия, начальные условия, сеточно-характеристические расчетные схемы. Затем на их основе были написаны модуль к языку Python, в результате чего появилась возможность работать с классами из языка С++ прямо из языка Python. Такой подход сильно упрощает отладку приложений, поскольку запуск можно вести в интерактивном режиме, используя оболочку языка Pyton. Помимо этого часть кода, отвечающая за организацию структуры расчетного модуля, частичную подготовку начальных данных, задание геометрии и структуры блоков, частично сохранение вынесены с систему скриптов.

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

Препроцессинг

Препроцессинг включает в себя следующие этапы:

• задание расчетной геометрии;

• задание граничных условий;

• задание параметров области начального возмущения;

• задание условий на границе раздела сред и на границах соседних блоков;

• задание механических характеристик среды;

• задание вычислительного алгоритма;

• задание параметров сохранения.

Рассмотрим данные шаги подробнее.

Расчетная геометрия может задаваться несколькими способами. Прежде всего можно прямо в конфигурационном файле указать параметры всех расчетных сеток, а именно число узлов по всем осям, шаги сеток по всем осям и их положение в пространстве. Криволинейные сетки могут быть загружены из файлов, на данный момент поддерживается загрузка таких сеток их текстового варианта формата УТК, однако за счет использования скриптового языка, в программу легко добавить новый функционал по загрузке сеток, и

Глава 5

Результаты численного моделирования ряда задач

Волны Лява Описание

В слоистых средах возможно возникновение определенных типов волн -волн Лява. Вектор смещения у таких воли параллелен границе раздела сред и перпендикулярен направлению распространения, т. е. волны Лява имеют горизонтальную поляризацию. В отличии от волн Рэлея, возникающих в одном полупространстве со свободной границей, волны Лява возникают в структурах типа упругий слой на упругом полупространстве. Теорию этих волн дал Ляв в 1911 г., именем которого они и названы.

Получить выражение для скорости волны в явном виде проблематично [31], определим волну по косвенным признакам и проверим, что соотношение скорости и длинны волны удовлетворяют уравнению волны Лява. Обозначим через cl скорость распространения волны Лява. Тогда для нее выполнено соотношение:

Cil <CL< Cf2, (5.1) где Cii и Ct2 - поперечные скорости звука для верхнего слоя и полупространства соответственно. Отсюда в частности следует, что для существования волн Лява необходимо, чтобы поперечная скорость звука в веществе слоя была меньше поперечной скорости звука в веществе полупространства. Также, в отличии от волн Рэлея, волны Лява обладают дисперсией, т. е. зависят от частоты и не сохраняют форму импульса. Амплитуда волны в полупростран

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

Волны Рэлея Описание

Волны Рэлея возникают на границе полуплоскости, заполненной однородной изотропной упругой средой. Теоретически эти волны были открыты Рэлеем в 1885 году, и они могут существовать вблизи свободной границы твердого тела, граничащей с вакуумом или достаточно разреженной газовой средой. Фазовая скорость таких волн направленна параллельно границе, а колеблющиеся вблизи нее частицы среды имеют как поперечную, перпендикулярную поверхности, так и продольную составляющие вектора смещения. Фазовая скорость волн Рэлея не зависит от длинны волны, т. е. они являются без дисперсными. Эти волны очень быстро затухают по глубине полуплоскости [32], вследствие наличия экспоненциальных множителей с показателем —где д - волновое число, г - координата, направленная вглубь полуплоскости, а - некий множитель, зависящий от параметров среды и скорости волны Рэлея. Отсюда следует, что чем меньше длина волны (больше волновое число), тем быстрее происходит затухание. Получается, что волны Рэлея поверхностные, т. е. их основная энергия сосредоточена у границы.

Волны Рэлея имеют большое значение в сейсмологии, поскольку наблюдаются вдали от эпицентра землетрясений, и именно они являются основной причиной разрушения наземных объектов.

Обозначим через ср, с^ - продольную и поперечную скорость в среде соответственно, а через сд - скорость волны Рэлея. Соотношение между этими скоростями задается уравнением третей степени [31]: = С3 - 8£2 + 8(2 + х)£ - 8(1 + X) = о, (5.5) где

Как известно, для всех упругих тел справедливо неравенство 0 < г) < \/2, при учете этого анализ (5.5) показывает, что 0,8741 < < 0, 9554. Таким образом, скорость волн Рэлея мало отличается от скорости сдвиговых волн, но всегда меньше ее.

Можно выписать явный вид корней уравнения (5.5) при некоторых частных значениях упругих постоянных среды [32]:

1. х = О, V2 = 2, <£ - 3 - ч/б.

2. Х = 1/3, V2 = 3, £ = 2(1 - - среда Коши.

Результаты моделирования

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

Расчетная область для всех тестов представляла собой параллелепипед размерами 1500 х 500 х 200 м, соответственно по осям х, у иг. На верхней грани области по оси г было выставлено граничное условие свободной границы, на других гранях устанавливалось граничное условие поглощения. Шаг сетки во всех расчетах брался порядка 10 м, число узлов в сетке - около 160 тыс. Шаг интегрирования по времени был выбран исходя из выполнения условия Куранта, во всех тестах он равен 0,001615 с. На рис. 5.5 представлена

Рис. 5.13. Расчетная геометрия скорость и доходит до земной поверхности позднее.

На рис. 5.15 показаны сейсмограммы, взятые на дневной поверхности у данных расчетов. Хорошо видно, что от взрыва приходит сильный фронт продольной волны, а от землетрясений сильный фронт поперечной. Полученный результат имеет хорошее сходство с экспериментальными данными. Также на сейсмограмме хорошо заметны отражения от всех слоев земной породы.

Прохождение волн через водонасыщенные грунты Описание

В данном разделе рассмотрена задача прохождения сейсмических воли через водонасыщенные грунты. Работа выполнялась совместно со специалистами из Российского Федерального Ядерного Центра — Всероссийского научно-исследовательского института экспериментальной физики (ФГУП "РФЯЦ-ВНИИЭФ").

Специалистами из Сарова был создан программный комплекс "Нимфа"для численного моделирования задач фильтрации. Этот комплекс ориентирован на решение традиционных гидрогеологических задач, связанных с прогнозом техногенного режима подземных вод при работе водозаборов и ведении в виде фронта упругой продольной волны, направление фронта - к земной поверхности.

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

1.0173 0Я16

Рис. 5.18. Распределение модуля скорости на поверхности земли в моменты времени i = 0.131,0.138,0.146,0.153,0.160,0.167 с от начала расчета

На рис. 5.19 показано распределение модуля скорости во вмещающем массиве.

Как видно из расчетов, различная водонасыщенность вносит существенные неоднородности в картинку распространения упругих волн. На дневной поверхности земли возникли сейсмические поверхностные волны, их распространение хорошо видно по рис. 5.18.

Проведена сшивка результатов расчета программного комплекса "Нимфа" для численного моделирования задач фильтрации, с программным комплексом, разработанным в рамках данной работы. На тестовых расчетах показано влияние насыщенности пор водой на изменение в распространении упругих волн.

1. Application of MPI-technology for solving the Boltzmann equation / N. I. Khokhlov, Yu. Yu. Kloss, B. A. Shurigin, F. G. Tcheremissine // 20-th International Conference on Transport Theory. — Obninsk: 2007. — P. 189.

2. Backus G. E. Long-Wave Elastic Anisotropy Produced by Horizontal Layering // Journal of Geophysical Research. — 1962. — oct. — Vol. 67. — Pp. 4427-4440.

3. Bermudez A., Hervella-Nieto L., Rodriguez R. Finite element computation of three-dimensional elastoacoustic vibrations // Journal of sound and vibration. 1999. - Vol. 219, no. 2. - Pp. 279-306.

4. Bohlen T., Saenger E.H. Accuracy of heterogeneous staggered-grid finite-difference modeling of Rayleigh waves // Geophysics. — 2006. — Vol. 71. — P. T109.

5. Chakravarthy S.R., Osher S. High resolution applications of the Osher upwind scheme for the Euler equations // 6th Computational Fluid Dynamics Conference. 1983. - Pp. 363-372.

6. Courant R.} Isaacson E., Rees M. On the solution of nonlinear hyperbolic differential equations by finite differences // Communications on Pure and Applied Mathematics. — 1952. — Vol. 5, no. 3. — Pp. 243-255.

7. Cronin VS. A draft primer on focal mechanism solutions for geologists: On the Cutting Edge Workshop on Teaching Structural Geology in the 21st Century // Science Education Resource Center, Carleton College. — 2004.

8. Dumbser M., Käser M. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes-II. The three-dimension

9. Lien F. S., Leschziner M. A. Upstream monotonic interpolation for scalar transport with application to complex turbulent flows // International Journal for Numerical Methods in Fluids. — 1994. — sep. — Vol. 19. — Pp. 527-548.

10. Nature of the scattered seismic response from zones of random clusters of cavities and fractures in a massive rock / VB Leviant, IB Petrov, FB Chel-nokov, IY Antonova // Geophysical prospecting. — 2007. — Vol. 55, no. 4.- Pp. 507-524.

11. Roe PL. Characteristic-based schemes for the Euler equations // Annual review of fluid mechanics. — 1986. — Vol. 18, no. 1. — Pp. 337-365.

12. Sibson R. A brief description of natural neighbour interpolation // Interpreting multivariate data. — 1981. — Vol. 21.

13. Sweby P. K. High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws // SIAM Journal on Numerical Analysis. — 1984. — oct.- Vol. 21. Pp. 995-1011.

14. Uniformly high order accurate essentially non-oscillatory schemes, III / A. Harten, B. Engquist, S. Osher, S.R. Chakravarthy // Journal of Computational Physics. — 1987. — Vol. 71, no. 2. — Pp. 231-303.

15. Virieux J. P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method // Geophysics. — 1986. — Vol. 51, no. 4. — Pp. 889-901.

16. Waterson NP, Deconinck H. A unified approach to the design and application of bounded higher-order convection schemes // Numerical methods in laminar and turbulent flow. — 1995. — Vol. 9. — Pp. 203-214.

17. Белоцерковский О. M., Агапов П. И., Петров И. Б. Численное моделирование последствий механического воздействия на мозг человека при

18. Паралелльные алгоритмы численных схем решения уравнения Больцма-на на основе технологии MPI / Е. П. Дербакова, Ю. Ю. Клосс, Н. И. Хохлов и др. // Электронный журнал "Исследовано в России". — 2007. — № 54. С. 581-588.

19. Петров И. Б., Тормасов А. Г., Холодов А. С. О численном изучении нестационарных процессов в деформируемых средах многослойной структуры // Известия АН СССР. Механика твердого тела. — 1989. № 4. - С. 89-95.

20. Петров И. Б., Холодов А. С. Численное исследование некоторых динамических задач механики деформируемого твердого тела сеточно-харак-теристическим методом // Журнал вычислительной математики и математической физики. — 1984. — Т. 24, № 5. — С. 722-739.

21. Петров И. Б., Хохлов Н. И. Моделирование сейсмических явлений се-точно-характеристическим методом // Труды МФТИ. — 2011. — Т. 3, № 3. С. 159-167.

22. Петров И. Б., Челноков Ф. Б. Численное исследование волновых процессов и процессов разрушения в многослойных преградах // Журнал вычислительной математики и математической физики. — 2003. — Т. 43, № 10. С. 1562-1579.

23. Программно-моделирующая среда для исследования течений газа в микро- и наноструктурах на основе решения уравнения Больцмана /

24. Ю. Ю. Клосс, Ф. Г. Черемисин, Хохлов Н. И., Шурыгин Б. А. // Атомная энергия. 2008. — Т. 105, № 4. - С. 211-213.

25. Разработка численных схем решения кинетического уравнения в кластерных средах на основе технологии MPI / Ю. Ю. Клосс, Ф. Г. Черемисин, Н. И. Хохлов, Б. А. Шурыгин // Информационные процессы. — 2007. — Vol. 7, по. 4. — Pp. 425-431.

26. Субботина А. Ю., Хохлов Н. И. Реализация клеточных автоматов "игра "Жизнь и клеточного автомата Кохомото-Ооно с применением технологии MPI // Компьютерные исследования и моделирование. — 2010. Vol. 2, по. 3. - Pp. 319-322.

27. Федоренко Р. П. Применение разностных схем высокой точности для численного решения гиперболических уравнений // Журнал вычислительной математики и математической физики. — 1962. — Vol. 2, по. 6. Pp. 1122-1128.

28. Холодов А. С. О построении разностных схем с положительной аппроксимацией для уравнений гиперболического типа // Журнал вычислительной математики и математической физики. — 1978. — Т. 18, № 6. — С. 1346-1358.

29. Холодов А. С. О построении разностных схем повышенного порядка точности для уравнений гиперболического типа // Журнал вычислительной математики и математической физики. — 1980. — Т. 20, № 6. — С. 1601-1620.

30. Хохлов Н. И. Тестирование латентности и скорости обмена данными библиотеки MPI при использовании сети Myrinet // Модели и методы обработки информации. Сборник научных трудов. — Москва: МФТИ, 2009.- С. 107-115.

31. Хохлов Н. И. Распараллеливание метода решения задач моделирования волновых процессов в деформируемом твердом теле используя технологию MPI // Математические модели и задачи управления. Сборник научных трудов. — Москва: МФТИ, 2011. С. 112-116.

32. Челноков Ф. Б. Численное моделирование деформационных динамических процессов в средах со сложной структурой: Ph.D. thesis / МФТИ (ГУ). 2005.

33. Чудов Л. А. Некоторые применения разностных методов в механике жидкостей и газа // Дисс. доктора физ.- матем. наук. М.: Ин- т проблем механ. АН СССР. — 1967.

34. Чушкин П.И. Метод характеристик для пространственных сверхзвуковых течений. Труды ВЦ АН СССР, 1968. - С. 121.

35. Шевченко A.A. Сейсмические исследования в скважинах. — МГУ. Геологический факультет. Кафедра сейсмометрии и геоакустики, 2007. Р. 136.