CC BY
102
УДК 519.633
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ УПРУГИХ ВОЗМУЩЕНИЙ, РАСПРОСТРАНЯЮЩИХСЯ ИЗ ОЧАГА ЗЕМЛЕТРЯСЕНИЯ1
В.И. Голубев, Н.И. Хохлов
Целью данной работы является разработка программного комплекса для проведения математического моделирования распространения упругих возмущений, зарождающихся в гипоцентре землетрясения, в гетерогенных средах. В качестве источника возмущения используется геофизическая модель очага землетрясения в виде подвижки по уже существующему разлому. Для описания динамического поведения среды используется система уравнений (гиперболического типа) механики твердого тела с явным выделением контактных границ неоднородностей. Она решается с помощью сеточно-характеристического метода на криволинейных структурных трехмерных сетках. Одной из особенностей используемого вычислительного алгоритма является его высокая масштабируемость по количеству вычислителей. Использование криволинейных сеток позволяет с высокой точностью описать большой класс расчетных геометрий. Авторами была проведена математическая постановка задачи, разработан программный комплекс и выполнена серия численных расчетов. В статье приведены результаты расчета распространения сейсмического возмущения в многослойной геологической среде и оценки сейсмостойкости наземного сооружения. Проведено исследование ускорения процесса расчета при увеличении количества вычислительных узлов.
Ключевые слова: сейсмическая активность, распространение упругих волн, гетерогенные среды, математическое моделирование, сеточно-характеристический численный метод, гиперболические системы уравнений.
Введение
Во всем мире ежегодно происходит более десяти крупных землетрясений, которые приводят к огромному числу человеческих жертв, множественным разрушениям строений и существенному материальному ущербу. Таким образом, оценка безопасности жилых и промышленных сооружений при воздействии на них динамических нагрузок является одной из актуальных проблем, как при выборе мест для строительства объектов на этапе проектирования, так и при последующем сопровождении конструкций с целью предупреждения чрезвычайных ситуаций природного характера. Под динамическими нагрузками здесь понимаются, как природные, так и техногенные землетрясения.
В настоящее время наибольшее распространение получила методика оценки сейсмостойкости, называемая нелинейный динамический анализ. Она основана на замене реальной трехмерной геометрии строения некоторой осредненной аналитической моделью. Для нее проводится оценка отклика на сейсмическое воздействие. Стоит отметить, что входными данными модели, как правило, являются экспериментальные акселерограммы (значения ускорения на поверхности Земли). Так, например, в работах [1, 2] проводится оценка сейсмостойкости мостов, в работах [3, 4] - жилых сооружений. Необходимо от-
1 Статья рекомендована к публикации программным комитетом Суперкомпьютерного форума «Суперкомпьютерные технологии в образовании, науке и промышленности 2012».
метить, однако, что точность таких оценок невысока в силу значительного упрощения физической модели процесса.
В настоящей работе разрабатывает подход, основанный на полном моделировании динамических процессов, происходящих в геологической среде при инициации сейсмической активности. Использование определяющей системы уравнений механики твердого тела совместно с реологическими соотношениями и явным выделением границ всех неоднородностей позволяет провести математическое моделирование с учетом реальной геометрии и свойств, как геологического массива, так и наземного строения. Для идентификации мест инициации разрушений возможно использование критериев наступления пластических или хрупких режимов деформирования.
Поскольку определяющая система уравнений имеет гиперболический тип, для ее численного решения используется сеточно-характеристический метод [5] на криволинейных структурных сетках. Данный метод хорошо себя зарекомендовал для решения задач сейсморазведки в существенно гетерогенных средах [6]. Также в работе [7] продемонстрирована возможность точного моделирования поверхностных волн (Релея и Лава), а в работе [8] - перспективы по применению данного метода для задач сейсмостойкости (в двумерной постановке).
В первом разделе подробно описана предложенная численная механикоматематическая модель очага землетрясения, основанная на геофизическом представлении о его природе. Перечислены параметры модели и предложен способ их восстановления по данным наземных наблюдений. Во втором разделе описана математическая модель геологической среды и численный метод, используемый в работе. В третьем разделе описан разработанный авторами программный комплекс и приведены результаты экспериментальной оценки масштабируемости реализованных вычислительных алгоритмов. В четвертом разделе приведены результаты моделирования землетрясения в слоистой гетерогенной среде, а также расчета сейсмостойкости наземного сооружения. В заключении приведено обобщение результатов работы.
1. Геофизическая модель очага землетрясения
Одной из особенностей данной работы является использование в качестве источника упругих возмущений механико-математической модели очага землетрясения. Отметим, что в геофизике существует множество различных моделей гипоцентра землетрясения [9], каждая из которых имеет экспериментальные факты, как свидетельствующие в пользу неё, так и против нее. Как правило, при возрастании сложности модели увеличивается число ее свободных параметров. Это приводит к усложнению однозначного решения задачи по восстановлению значений параметров модели по данным наземных измерений.
В настоящей работе было решено взять за основу максимально простую, но в то же время физически корректную модель, называемую «подвижка по разлому». Основой данной модели является геофизическое представление о землетрясении, как о процессе постепенного накопления упругих напряжений в земной коре, состоящей из отдельных блоков. Поскольку разломы между блоками заполнены брекчией, эффективная сила трения между бортами достаточно велика. Однако, при достижении некоторого порога, происходит проскальзывание граничащих блоков вдоль разлома и, соответственно, резкий спад напряжений. При этом вблизи разлома могут наблюдаться, как упругие, так и
пластические и хрупко-пластические процессы. На большом удалении от места инициации возмущения существенны лишь упругие процессы, которые выражаются в распространении упругих волн через геологическую среду.
Рассмотрим подробнее предлагаемую модель очага землетрясения. На рис. 1 изображен контакт двух блоков (линия контакта - разлом). Для задания ориентации плоскости разлома принято использовать два угла, определяемых следующим образом. Один из них называется strike angle (ft) и равен углу, который образуют линия пересечения плоскости разлома с горизонтальной плоскостью и направление на север. Второй угол, называемый dip angle (Ô), образуется между плоскостью разлома и горизонтальной плоскостью. Вдоль разлома происходит проскальзывание граничащих блоков, причем делается предположение о чистом продольном сдвиге, т.е. отсутствии составляющей скорости движения, перпендикулярной плоскости разлома. Ориентация вектора скорости задаётся углом между направлением вдоль линии пересечения плоскости разлома с горизонтальной плоскостью (strike direction) и направлением проскальзывания блоков (slip direction). Данный угол называется rake angle (Л). Кроме того, для задания интенсивности возмущения используется также абсолютное значение скорости смещения. При этом предполагается, что в некоторой области, прилегающей к разлому, модуль вектора скорости постоянен по пространству. При переходе через плоскость разлома направление скорости меняется строго на противоположное направление.
По результатам измерений сейсмометров, расположенных на дневной поверхности, возможно определение ориентации плоскости разлома. Для оценки величины начальной скорости в источнике возмущения может быть использовано сравнение результатов численного моделирования с амплитудами зарегистрированных колебаний.
Рис. 1. Модель очага землетрясения
2. Математическая модель среды и численный метод
Для математического моделирования волновых процессов в деформируемом твердом теле использовалась система динамических уравнений, объединяющая уравнения движения и реологические соотношения в виде:
Р ■ V =У J ■ °у, ау = Чщ ■ £к1 + Ру ■
Здесь Р - плотность среды, Vi - компоненты скорости смещения, а у , skl - компоненты тензоров напряжения и деформаций, V j - ковариантная производная по j-й координате, F у - добавочная правая часть.
Вид компонент тензора 4-го порядка qijkl определяется реологией среды. Для линейно-упругого тела
Qijkl = ^àij àkl + M (àik àjl + àil àjk )
Вторая группа уравнений представляет собой продифференцированный по времени закон Гука:
^ij = ^(& 11 + &22 + &33 ")^ij + 2M&ij
В этом соотношении X и ^ - упругие постоянные Ляме, 5у - символ Кронекера.
Для численного решения динамических уравнений механики деформируемого твердого тела использовался сеточно-характеристический метод, подробнее с которым можно ознакомиться в [5, 6]. Возможность описания произвольных геометрий реализована за счет использования структурных криволинейных сеток. При этом выделение неоднородностей (например, различных геологических слоев) производится уже на этапе построения расчетной сетки. При этом для сшивки решений на границе областей используется постановка контактных условий в явном виде, не приводящая к появлению паразитных переотраженных волн.
3. Описание расчетного программного комплекса
На основе описанных моделей и методов был разработан программный комплекс, позволяющий проводить компьютерное моделирование динамических процессов, происходящих в геологической среде и материалах строений, при инициации сейсмической активности. В данный комплекс входят модули препроцессинга, расчета и постпроцессинга.
Для реализации совместимости с open source проектами использованы форматы GMSH и RED представления структурных сеток. Таким образом, задание произвольной геометрии возможно с использованием GUI GMSH.
Расчетный модуль использует технологию MPI для использования большого числа вычислителей. Начальная расчетная сетка делится на возможно равные части между процессами. При этом каждый процесс хранит только свою часть сетки, что позволяет существенно уменьшить потребление памяти и проводить расчеты на сетках, содержащих миллиарды узлов. Каждый процесс хранит еще слой приграничных ячеек для обмена актуальными данными с соседними процессам при расчете значений на следующем временном слое. Размер приграничного слоя зависит от порядка используемой разностной схемы. Так, для второго порядка точности он должен содержать две ячейки. Эффективность такого распараллеливания была протестирована на серии расчетов с сеткой, содержащей 64 миллиона узлов. При этом сравнивалось время расчета на разном количестве вычислителей: от 1 до 108 включительно. На рис. 2 приведены полученные зависимости, свидетельствующие о том, что масштабируемость составляет порядка 80 %.
Оценка сейсмостойкости строений основана на идентификации мест инициации разрушений при действии динамической нагрузки. На каждом временном шаге произво-
дится перебор всех узлов сеток, описывающих наземное строение. В каждом узле проверяется выполнение критерия Мизеса пластических деформаций. Те точки, в которых критерий выполнен, помечаются маркером, свидетельствующим о начавшихся в них разрушениях. Таким образом, по результатам расчета строится картина зависимости зон разрушений в строении от времени.
Для возможности детального анализа результатов расчета реализован функционал сохранения полного тензора напряжений, маркеров разрушений и вектора скорости в формате УТК.
Рис. 2. Ускорение вычислительного алгоритма при увеличении числа потоков
4. Результаты численного моделирования
Одним из преимуществ компьютерного моделирования сейсмической активности является возможность получения полной волновой картины в виде зависимости тензора упругости и вектора скорости от времени. При этом, конечно, сохраняется возможность построения всех общепринятых форматов (сейсмограмм, отдельных сейсмотрасс) представления процесса для сравнения с результатами полевых измерений.
В настоящей работе было проведено численное моделирование распространения упругого возмущения из очага землетрясения через слоистую геологическую среду. Она состояла из четырех слоев с различными упругими характеристиками, значения которых приведены в таблице. Результаты численного моделирования приведены на рис. 3 (изображено вертикальное сечение) в виде последовательности волновых картин. Детальный анализ волновой картины [8] показывает, что из гипоцентра распространяются четыре волны - две продольные и две поперечные. Когда они достигают границы разде-
ла слоев, в следующий слой проходит «проходящая» волна, а в предыдущий слой возвращается «отраженная» волна. Таким образом, в общем случае происходит ослабление всего возмущения.
Таблица
Характеристики слоистой среды
Номер слоя Толщина слоя, м Плотность, кг/м3 Продольная скорость, м/с Поперечная скорость, м/с
1 300 2500 4190 2793
2 400 2500 4650 3100
3 500 2500 5250 3500
4 1600 2500 5850 3900
г
Рис. 3. Волновые картины распространения сейсмического возмущения
в многослойной среде
В настоящей работе также было проведено полное моделирование динамических процессов в области, включающей, как геологическую среду, так и наземное сооружение. Характерные размеры строения: 6 м * 6 м * 6 м. Толщина стен была выбрана равной 0,5 м. Рассматривалось влияние одной компоненты возмущения, распространяющегося из очага землетрясения, - поперечной волны большой амплитуды. Параметры вмещающего массива: плотность - 2000 кг/м , скорость распространения продольной волны - 2000 м/с, скорость распространения поперечной волны - 1200 м/с. Материал строения имеет следующие свойства: плотность - 2500 кг/м , скорость распространения продольной волны - 4000 м/с, скорость распространения поперечной волны - 2500 м/с.
Результатами моделирования являются рассчитанные поля скорости и напряжений, эволюционирующие во времени. Также был проведен поиск мест, в которых начинается разрушение здания. Расчетная сетка содержала порядка 13 миллионов узлов, а наибольший размер ячейки составлял 6 см. Результаты моделирования приведены на рис. 4.
Рис. 4. Модуль скорости (слева) и места разрушений (справа)
Заключение
В настоящей работе предложена численная механико-математическая модель гипоцентра землетрясения. Проведена разработка программного комплекса, позволяющего проводить полноволновое моделирование динамических процессов, происходящих в геологической среде при инициации сейсмической активности. Эмпирическая оценка масштабируемости свидетельствует о высокой степени параллелизма используемых программных алгоритмов. На основе серии расчетов изучено распространение упругих возмущений из очага землетрясения в слоистой геологической среде. На основе анализа выполнения критерия Мизеса пластических деформаций проведена идентификация мест разрушений в наземном строении путем прямого численного моделирования.
Одной из особенностей геологических сред является их внутренняя трещиноватость: наличие газонасыщенных и флюидонасыщенных трещин протяженностью от нескольких сантиметров до десятков метров. Для корректного моделирования волновых процессов вблизи таких неоднородностей необходимо использовать очень подробную сетку. Однако равномерное измельчение сетки во всей расчетной области существенно увеличит время расчета. Таким образом, дальнейшим направлением исследования предполагается разработка алгоритмов сшивки структурных и неструктурных сеток, позволяющих описать неоднородности с использованием локального измельчения расчетной сетки.
Исследование выполнено при поддержке Министерства образования и науки Российской Федерации, соглашение 14-132.21.1809, и РФФИ в рамках научного проекта № 12-07-31028.
Литература
1. Ala Saadeghvaziri, M. Seismic behavior and capacity/demand analysis of three multispan simply supported bridges / M. Ala Saadeghvaziri, A.R. Yazdani-Motlag / / Engineering Structures. - 2008. - Vol. 30. - P. 54-66.
2. Soneji, B.B. Influence of soil structure interaction on the response of seismically isolated cable-stayed bridge / B.B. Soneji, R.S. Jangid // Soil Dynamics and Earthquake Engineering. - 2008. - Vol. 28. - P. 245-257.
3. Maniyar, M.M. Probabilistic seismic performance evaluation of non-seismic RC frame buildings / M.M. Maniyar, R.K. Khare, R.P. Dhakal / / Structural Engineering and Mechanics. - 2009. - Vol.33, No. 6. - P. 725-745.
4. Solberg, K.M. Rapid expected annual loss estimation methodology for structures / K.M. Solberg, R.P. Dhakal, J.B. Mander, B.A. Bradley / / Earthquake Engineering and Structural Dynamics. - 2008. - Vol.37, No. 1. - P. 81-101.
5. Магомедов, К.М. Сеточно-характеристические численные методы / К.М. Магомедов, А.С. Холодов - М.: Наука, 1988. - 288 с.
6. Квасов, И.Е. Численное моделирование сейсмических откликов в многослойных геологических средах сеточно-характеристическим методом / И.Е. Квасов, С.А. Панкратов, И.Б. Петрова / / Математическое моделирование. - 2010. - Т. 22, № 9. -
С. 13-22.
7. Хохлов, Н.И. Моделирование сейсмических явлений сеточно-характеристическим методом / Н.И. Хохлов, И.Б. Петров // Труды МФТИ. - 2011. - Т. 3. - С. 159-167.
8. Голубев, В.И. Воздействие природных катастроф на наземные сооружения /
В.И. Голубев, И.Е. Квасов, И.Б. Петров // Математическое моделирование. - 2011.
- Т. 23, № 8. - С. 46-54.
9. Поле упругих напряжений Земли и механизм очагов землетрясений / Л.М. Балаки-на, А.В. Введенская, Н.В. Голубева и др. - М.: Наука, 1972. - 191 с.
Голубев Василий Иванович, аспирант, Московский физико-технический институт, w.golubev@ mail. ru
Хохлов Николай Игоревич, инженер, Московский физико-технический институт, k h@inbox.ru
MATHEMATICAL MODELLING OF ELASTIC PERTURBATIONS PROPAGATING FROM THE EARTHQUAKE HYPOCENTER
V.I. Golubev, Moscow Institute of Physics and Technology (Dolgoprudny, Russian Federation),
N.I. Khokhlov, Moscow Institute of Physics and Technology (Dolgoprudny, Russian Federation)
The goal of this article is the development of the research software for the mathematical modelling of elastic waves propagation originated in the earthquake hypocenter through heterogeneous media. As a source of the perturbation the geophysical focal mechanism model based on the slipping along the fault is used. For the description of dynamic behavior of media the hyperbolic system of equations of elastic media is used with the explicit allocation of contact borders of heterogeneities. It is solved using the grid-characteristic method on curvilinear structured 3D meshes. One feature of used numerical algorithm is its high scalability per core. The usage of curvilinear meshes allows describing a wide variety of geometries with high precision. Mathematical formulation of problem, development of the research software and a set of numerical experiments were done by authors. The results of modelling of propagation of seismic perturbation through geological multilayered medium and assessment of seismic resistivity of ground facility are described in this article. The estimation of software scalability per core was carried out.
Keywords: seismic activity, propagation of elastic waves, heterogeneous media, mathematical modeling, numerical grid-characteristic method, hyperbolic system of equations.
References
1. Ala Saadeghvaziri M., Yazdani-Motlag A.R. Seismic behavior and capacity/demand analysis of three multi-span simply supported bridges // Engineering Structures. 2008. Vol. 30. P. 54-66.
2. Soneji B.B., Jangid R.S. Influence of soil structure interaction on the response of seismi-cally isolated cable-stayed bridge // Soil Dynamics and Earthquake Engineering. 2008. Vol. 28. P. 245-257.
3. Maniyar M.M., Khare R.K., Dhakal R.P. Probabilistic seismic performance evaluation of non-seismic RC frame buildings // Structural Engineering and Mechanics. 2009. Vol.33, No. 6. P. 725-745.
4. Solberg K.M., Dhakal R.P., Mander J.B., Bradley B.A. Rapid expected annual loss estimation methodology for structures // Earthquake Engineering and Structural Dynamics. 2008. Vol.37, No. 1. P. 81-101.
5. Magomedov K.M., A.S. Holodov Setochno-haracteristicheskie chislennie metodi [Numerical Grid-Characteristic Methods]. M.: Science, 1988. 288 p.
6. Kvasov I.E., Pankratov S.A., Petrov I.B. Chislennoe modelirovanie seismicheskih ot-klikov v mnogosloinih geologicheskih sredah setochno-haracteristicheskim metodom [Numerical Modelling of Seismice Responses in Multilayered Geological Media using Grid-Characterisitic Methods]. Matematicheskoe modelirovanie [Mathematical Modelling]. 2010. Vol. 22, No 9. P. 13-22.
7. Khokhlov N.I., Petrov I.B. Modelirovanie seismicheskih yavlenij setochno-haracteristicheskim metodom [Modelling of Seismic Phenomena using Grid-Characteristic Method]. Trudi MFTI [Proceedings of MIPT]. 2011. V. 3. P. 159-167.
8. Golubev V.I., Kvasov I.E., Petrov I.B. Influence of Natural Disasters on Ground Facilities // Mathematical Modelling and Computer Simulations, 2012. Vol. 4, No. 2. P. 129134.
9. Balakina L.M., Vvedenskaya A.V., Golubeva N.V., Misharina L.A., Shirokova E.I. Pole uprugih napryazhenij Zemli i mechanismi ochagov zemletryasenij [Elastic Stress Field of the Earth and earthquake focal mechanisms]. M.: Science, 1972. 191 p.
Поступила в редакцию 5 марта 2013 г.
