Моделирование сейсмических явлений сеточно-характеристическим методом Текст научной статьи по специальности «Физика»

УДК 519.63

Н.И. Хохлов, И.Б. Петров Московский физико-технический институт (государственный университет)

Моделирование сейсмических явлений сеточно-характеристическим методом

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

Ключевые слова: сеточно-характеристический метод, волны Рэлея, волны Лява, динамика деформируемого твердого тела, моделирования, системы гиперболических уравнений.

I. Введение

Математическое моделирование сейсмостойкости и сейсмических явлений имеет важное значение на этапе проектирования и выбора площадок для постройки зданий, мостов, плотин, атомных электростанций, других стратегических объектов и сложных наземных сооружений. В данной работе рассмотрен подход моделирования таких явлений используя сеточно-характеристический метод решения уравнений упругости на параллелепипедных сетках в трехмерном случае. Важными свойствами данного метода является монотонность, возможность строить корректные численные алгоритмы на контактных границах и границах области интегрирования. При численном решении весьма специфических задач геодинамики, имеющих ярко выраженный волновой характер с большим количеством контактных и свободных границ, данный метод позволяет получать адекватные численные решения, моделирующие сложные динамические процессы, происходящие в существенно неоднородных средах. В [1, 2] этот метод исследовался для численного моделирования динамических процессов в упругопластических средах на прямоугольных сетках в двумерном случае, в работах [3, 4] — на треугольных.

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

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

II. Математическая модель

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

Сформулируем основные уравнения линейной динамической теории упругости, которым подчиняется состояние бесконечно малого объема линейно-упругой среды. Они объединяют в себе уравнения движения и реологические соотношения в виде

ру = ^,

Т = А(У ■ V)I + р(У ® V + V ® V),

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

Рассмотрим приведенные нестационарные уравнения теории упругости для случая трех переменных к некоторой ортонормированной системе координат (жі,ж2,жз). Первая строка в системе уравнений (1) дает три уравнения движения, а вторая — шесть реологических соотношений. Составим вектор искомых функций, состоящий из 9-ти компонент:

и = {г>1 ,г2,г3,о11,о12,о13,о22,о23,о33}Т,

здесь гі — компоненты вектора скорости смещения, 07 — компоненты симметричного тензора напряжений Т. Тогда перечисленные модели твердого тела допускают запись системы уравнений (1) динамики деформируемого твердого тела в матричном виде [5]:

ди=£ а £ • (2)

.7 = 1 •’

где А7 — матрицы размера 9 х 9.

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

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

Для численного моделирования задач динамики деформируемого твердого тела широко применяется сеточно-характеристический метод [6, 7]. Вначале применяется метод расщепления по пространственным координатам, в результате чего имеем три одномерных системы:

ди 4 ди . .

« = А «7 0 =1А3- (3)

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

«и о-1 л о «и

ді = П ЛЛ дх3 •

где матрица П — матрица, составленная из собственных векторов, Лj — диагональная матрица, элементами которой являются собственные значения. Для всех координат матрица Л имеет вид

Л = б,гад{с1, — с1,с2, — с2,с2, — с2, 0, 0, 0},

где с1 = л/(Л + 2ц)/р — продольная скорость звука в среде, с2 = \/Д7Р — поперечная скорость звука.

После замены переменных V = Пи каждая из систем (3) распадается на девять независимых скалярных уравнений переноса (индекс о далее опускается. где это возможно):

«V . дv

т + Л дХ = °'

Одномерные уравнения переноса решаются с помощью метода характеристик либо обычными конечно-разностными схемами. Конкретно в данной работе все расчеты проводились используя ТУБ-схемы второго порядка точности с ограничителем вирегЪее [8].

После того, как все компоненты V перенесены, восстанавливается само решение:

ип+1 = П-1vп+1.

Н.Э. Расчетная сетка

Рассмотрим расчетную область ^ С Я3, задача состоит в восстановлении полей скорости и тензора напряжения на области ^ по заданному начальному распределению и граничным условиям. Структура среды известна и задается набором из т параллелепипедных блоков с ребрами, параллельными осям декартовой системы координат: ^ = ит=1 . Любой из блоков ^ характе-

ризуется своим набором параметров среды, граничных и контактных условий. В каждом блоке строятся прямоугольные сетки с постоянным внутри блока шагом = {к^1 ,^&2,к^3}, при этом в соседних блоках при условии наличия контактной границы сетки должны быть согласованы. Под согласованием сеток понимается равенство шагов сетки в плоскости контакта и совпадение координат узлов в прилегающих сетках. На рис. 1 изображены примеры расчетных областей. Для расчета на границах сеток в соседних блоках реализованы различные виды контактных границ.

Рис. 1. Пример расчетных областей

III. Волны Рэлея

Описание

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

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

Обозначим через ср,с^ — продольную и поперечную скорость в среде соответственно, а через с# — скорость волны Рэлея. Соотношение между этими скоростями задается уравнением третьей степени [10]:

/ (С) = С3 — 8£2 + 8(2 + х)С — 8(1 + х) = 0, (4)

где

2 2 с

С = с#/с2 ,х = 1 — - П = ср.

' I с

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

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

1. х = 0, п2 = 2, С = 3 - л/5.

2. х = 1/3, п2 = 3, С = 2(1 — ) — среда Коши.

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

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

Расчетная область для всех тестов представляла собой параллелепипед размерами

1500 х 500 х 200 м соответственно по осям х, у и г. На верхней грани области по оси г было

выставлено граничное условие свободной границы, на других гранях устанавливалось граничное условие поглощения. Шаг сетки во всех расчетах брался порядка 10 м, число узлов в сетке — около 160 тыс. Шаг интегрирования по времени был выбран исходя из выполнения условия Куранта, во всех тестах он равен 0,001615 с. На рис. 2 представлены расчетная область и зона задания начального возмущения. Во всех областях изначально заданы нулевые распределения скорости и напряжения, за исключением небольшого возмущения скорости на границе области по оси х. Вектор скорости возмущения направлен по оси г, величина — 0,001 м/с.

Первый расчет был проведен со следующими параметрами среды: с* = 2000 м/с

и ср = л/2с* « 2830 м/с. Расчетная скорость волн Рэлея в данном случае будет

сд = \^3 — л/5с* ~ 0,874с* ^ 1748 м/с. На рис. 3 представлено сечение, проходящее через середину расчетной области, параллельно осям х и г и перпендикулярно оси у. Отображены модули скорости для моментов времени £ = 0,2955 и £ = 0,6412 с.

На картинках хорошо видно распространение поверхностной волны, видно, что с глубиной проникновения волна сильно затухает, максимальная амплитуда находится у поверхности среды. За это время волна распространилась на 600 м, нетрудно посчитать ее фазовую скорость: с = 1736 м/с. Эта величина хорошо соотносится с теоретическим значением, полученным ранее (1748 м/с).

Во втором расчете поперечная скорость была такой же, а продольная ср = ^/с* ~ 3464 м/с, теоретическая скорость волны Рэлея в таком случае будет сд = ^2(1 — -^=)с* ~ 0,9194с* ^ 1839 м/с. На рис. 4 представлены аналогичные сечения для моментов времени £ = 0,2794 и £ = 0,6089 с.

Волновая картина аналогична картине, полученной в первом случае, так же несложно получить значение скорости волны: с = 1821 м/с, что снова хорошо согласуется с теоретическим значением.

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

Рис. 3. Распространение волны Рэлея в первом Рис. 4. Распространение волны Рэлея во втором

случае: а) £ = 0,2955 с, Ъ) £ = 0,6412 с случае: а) £ = 0,2794 с, Ъ) £ = 0,6089 с

IV. Волны Лява

Описание

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

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

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

с*1 < сь < с*2, (5)

где с*1 и с*2 — поперечные скорости звука для верхнего слоя и полупространства соответственно. Отсюда, в частности, следует, что для существования волн Лява необходимо, чтобы поперечная скорость звука в веществе слоя была меньше поперечной скорости звука в веществе полупространства. Также, в отличие от волн Рэлея, волны Лява обладают дисперсией, то есть зависят от частоты и не сохраняет форму импульса. Амплитуда волны в полупространстве затухает по экспоненте.

Обозначим через Н толщину верхнего слоя, и — угловая частота волны Лява, х1 = и2/с! ,х2 = и2/с22. Тогда конечное число волн Лява (количество возможных гармоник) в поверхностном слое определяется соотношением

N =[Н /х1— х2/п] +1

(6)

здесь [а] означает целую часть числа а. Подставив в (6) выражение для Х1 и Х2, а также с учетом того, что и = 2псь/Аь, Аь — длина волны Лява, получим

N =

2Нсь 11 Аь

41

1

42

+ 1.

Выпишем уравнение, определяющее свойства поверхностных волн Лява:

^ п = (М2/М1 )[(Х1Н )2 — (Х2Н )2 — п2 ]1/2 п

(7)

(8)

здесь Д2, М1 — упругие постоянные Ламе для полупространства и слоя соответственно, п = \/Х1 — С2Н, С2 = и2/сЬ. Данное уравнение определяет соотношение между скоростью и длиной волны Лява.

*

а

•

Ъ

*

с

ш щ

а

волн Лява

Рис. 6. Распространение волны Лява: а) £ = 0,258 с, Ъ) £ = 0,347 с, с) £ = 0,436 с, а) £ = 0,525 с

2

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

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

Расчетная область представляла собой два параллелепипеда, представляющие верхний слой и нижнее полупространство, их размеры 1500 х 500 х 50 м (Н = 50 м) и 1500 х 500 х 150 м соответственно по осям х, у и г. На верхней грани слоя по оси г было выставлено граничное условие

свободной границы, на других гранях слоя и полупространства устанавливалось граничное условие поглощения. Между слоем и полупространством установлено условие контактной границы. Шаг сетки был равен 10 м, число узлов в сетке — около 170 тыс. Шаг интегрирования по времени был выбран исходя из выполнения условия Куранта, в данном тесте он равен 0,001615 с. На рис. 5 представлена расчетная область и зона задания начального возмущения. Во всей области изначально задавалась нулевая скорость и напряжения, а на границе верхнего слоя по оси х задается небольшое возмущение скорости. Вектор скорости направлен по оси у, величина возмущения —

0,001 м/с.

Рис. 7. Компонента у скорости среды у поверхности слоя: а) £ = 0,258 с, Ь) £ = 0,347 с, с) £ = 0,436 с, -г, о - тт

_ кос Рис. 8. Определение угловой частоты волны Лява

О) Т — и,525 с

В расчете использовались следующие параметры среды (индекс 1 относится к верхнему слою, 2 — к полупространству): ср1 — 3000 м/с, вц — 2000 м/с, ср2 — 6000 м/с, с^2 — 3000 м/с, Р1 — р2 — 2500 кг/м3. На рис. 6 представлено сечение, проходящее через середину расчетной области, параллельно осям х и г и перпендикулярно оси у. Цветом отображены у-компоненты скорости для моментов времени 0,258, 0,347, 0,436 и 0,525 с.

Как видно из графиков, полученная волна обладает горизонтальной поляризацией, очень быстро затухает в полупространстве. Распространение происходит вдоль поверхностного слоя. По мере распространения форма импульса меняется. Несложно посчитать скорость полученной волны: с — 600/(0,525 — 0,258) — 2247 м/с, как видно, она удовлетворяет соотношению (5).

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

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

Подставив в формулу (7) полученные значения: скорость распространения, длину волны и параметры расчетной среды, получим число волн Лява: N — [0,5992] + 1 — 1, это говорит о том, что должна существовать только одна гармоника и уравнение (8) имеет только один действительный корень.

Уравнение (8) трансцендентное, и получить в явном виде значение для угловой частоты волны Лява невозможно. Найдем приближенное значение угловой частоты волны Лява из уравнения графически, предварительно подставив в него параметры упругих сред слоя и полупространства и полученное значение скорости волны Лява. Графики левой и правой частей уравнения (8) представлены на рис. 8.

По оси х отложено значение ш — угловая частота волны Лява. Данные графики пересекаются в одной точке и на рисунке представлена только часть графиков в соответствующем масштабе. То, что данное уравнение имеет один корень, хорошо соотносится с полученным ранее результатом. Как видно из графика, значение угловой частоты порядка еог — 107,7 Гц. Ранее уже было получено значение длины волны Лява А^ — 140 м, отсюда получим экспериментальное значение частоты: шехр — 2пс^/А^ — 101 Гц. Видно, что теоретическое значение отличается от экспериментального не более чем на 10%.

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

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

Описание

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

Рис. 9. Область расчета прохождения волн через здание

Расчетная область изображена на рис. 9 и представляет собой часть земли прямоугольной формы и размерами 20 х 20 х 10 м. На поверхности земли находится двухэтажное здание размерами порядка 5 х 5 х 6 м, фундамент углублен в земные породы на величину порядка 1 м. Параметры земли, при которых велся расчет: ср — 2000 м/с, с* — 1300 м/с, р — 2000 кг/м3. Здание рассчитывалось как монолитное, сделанное из бетона с параметрами среды: ср — 4000 м/с, с* — 2500 м/с, р — 4000 кг/м3, ау — 15 МПа. Здесь ау — предел текучести материала.

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

а- = л/э^ - у - ^22)2 + (^22 - аэ3)2 + (аэ3 - аі1)2 + 6(аі2 + а23 + азі)

здесь аі — интенсивность напряжений, а ^ — второй инвариант девиатора тензора напряжений. Исходя из условия Мизеса, определялось, в каких точках здания возможны разрушения. Бетон считался абсолютно не пластичным и разрушения происходили при выполнении условия а і > а у. Следует учесть, что производилось только определение узлов расчетной сетки, где выполнилось условие пластичности Мизеса, но никакие изменения в систему уравнений не вносились.

В области земли, везде кроме поверхности, были выставлены условия поглощения, на поверхности — свободная граница. На здании везде было условие свободной границы. На рис. 9 черным цветом показано начальное возмущение, в зависимости от расчета это продольная или поперечная волна. Модуль скорости возмущения во всех расчетах одинаковый и равен 0,9м/с.

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

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

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

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

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

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

Рис. 10. Результаты расчета для продольной вол- Рис. 11 Результаты расчета для поперечноИ в°л ны: а) Т — 0,0013 с, Ь) Т — 0,0039 с, с) Т — 0,0079 с ны: а) Т — 0,0020 с, Ь) Т — °,°°39 с, с) Т — °,°186 с

У.2. Заключение

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

Литература

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

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

3. Квасов И.Е., Панкратов С.А., Петров И.Б. Численное моделирование сейсмических откликов в многослойных геологических средах сеточно-характеристическим методом // Математическое моделирование. — 2010. — Т. 22, вып. 9. — С. 13-22.

4. Белоцерковский О.М., Агапов П.И., Петров И.Б. Численное моделирование последствий механического воздействия на мозг человека при черепно-мозговой травме // Журнал вычислительной математики и математической физики. — 2006. — Т. 46, вып. 3. — С. 1711-1720.

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

6. Магомедов К.М., Холодов А.С. О построении разностных схем для уравнений гиперболического типа на основе характеристических соотношений // Журнал вычислительной математики и математической физики. — 1969. — Т. 9, вып. 2. — С. 373-386.

7. Магомедов К.М., Холодов А.С. Сеточно-характеристические численные методы. — М.: Наука, 1988.

8. P.L. Roe. Characteristic-Based Schemes for the Euler Equations // Annual Review of Fluid Mechanics. — 1986. — N.18. — P. 337-365.

9. Горшков А.Г., Медведский А.Л., Рабинский Л.Н., Тарлаковский Д.В. Волны в сплошных средах. — М.: Физматлит, 2004.

10. Бреховских Л.М., Гончаров В.В. Введение в механику сплошных сред. — М.: Наука, 1982.

Поступила в редакцию 24.05.2011.