Журнал вычислительной математики и математической физики, 2020, T. 60, № 9, стр. 1546-1565

ВВЕДЕНИЕ

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

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

Задача рассмотрена при условии плоской деформации. Приведены определяющие гиперболические уравнения динамической теории упругости, которые описывают движение среды и штампа, краевые и контактные условия на границах и контактной поверхности. Для решения задачи принята явная разностная схема, построенная на основе метода бихарактеристик [1], [2] с привлечением идеи расщепления по пространственным координатам [3]–[5]. Дано описание точечной расчетной схемы и шаблона для изотропной упругой полуплоскости с прямоугольным упругим штампом на ее границе. Получены разрешающие разностные уравнения для внутренних, граничных, угловых, особых и контактных точек сопряжения штампа и полуплоскости.

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

1. ПОСТАНОВКА КОНТАКТНОЙ КРАЕВОЙ ЗАДАЧИ

Рассматривается полуплоскость ${{x}_{1}} \geqslant 0$ упругой однородной изотропной среды (первая среда ${{D}_{1}}$) с плотностью ${{\rho }^{{(1)}}}$ и упругими коэффициентами Ламе ${{\lambda }^{{(1)}}}$ и ${{\mu }^{{(1)}}}$, на поверхности которого расположено упругое тело (штамп – вторая среда ${{D}_{2}}$) высотой ${{d}_{1}}$ и шириной $2{{d}_{2}}$ с плотностью ${{\rho }^{{(2)}}}$ и коэффициентами Ламе ${{\lambda }^{{(2)}}}$, ${{\mu }^{{(2)}}}$.

До начального момента времени среды находились в статическом состоянии, которое здесь не рассматривается. Предполагается, что в начальный момент времени $t = 0$ в упругой полуплоскости на глубине $L$ возникает горизонтальная трещина $S$ (${{x}_{1}} = L$, $\left| {{{x}_{2}}} \right| \leqslant d$) (фиг. 1). При $t > 0$ динамическое состояние сред описывается векторами ${{{\mathbf{u}}}^{{(k)}}}$ ($k = 1,\;2$) (смещение $k$-й среды от статического состояния) и тензорами динамических напряжений $\sigma _{{ij}}^{{\left( k \right)}}$. Следовательно,

Границы полуплоскости (дневная поверхность) и штампа свободны от нагрузок: Условия на контактной границе отвечают требованиям жесткого контакта (полное сцепление): Так как на бесконечности отсутствуют источники колебания, то очевидным является требование, чтобы на бесконечности выполнялись условия затухания: При описанных условиях необходимо исследовать напряженно-деформированное состояние неоднородной среды ${{D}_{1}} \cap {{D}_{2}}$ при $t > 0$.

2. УРАВНЕНИЯ ДВИЖЕНИЯ УПРУГИХ ПОЛУПЛОСКОСТИ И ШТАМПА

Для решения поставленной задачи наряду с начальными (1.1) и граничными условиями (1.2)–(1.5) используется система дифференциальных уравнений, состоящая из уравнений движения упругой среды и закона Гука [6]:

Здесь $u_{i}^{{\left( k \right)}}$ – компоненты вектора перемещений, $\sigma _{{ij}}^{{\left( k \right)}}$ – тензор напряжений, $F_{i}^{{\left( k \right)}}$ – компоненты объемной силы, ${{\delta }_{{ij}}}$ – символ Кронекера.

Компоненты ${{F}_{i}}$ объемной силы $F$ определяются сингулярной обобщенной функцией – простым слоем на трещине $S$ [7]. Здесь они имеют следующий вид:

где выражение в квадратных скобках – скачок компонент тензора напряжений на берегах трещины, ${\mathbf{n}}$ – единичная нормаль к ее поверхности (в данном случае ${\mathbf{n}} = ({{n}_{1}},{{n}_{2}}) = (1;\;0)$), $H(x)$ – функция Хевисайда, $\delta (x)$ – функция Дирака. Предполагается, что скачок напряжений на трещине известен: Здесь в расчетах рассматривается импульс вида А для представления $\delta $-функции Дирака используются дельтообразные последовательности ${{\delta }_{\varepsilon }}(x)$ [8] где $\varepsilon > 0$, т.е. в расчетах

Для решения поставленной задачи используется метод бихарактеристик [1]–[5], на основе которого разработаны алгоритм и пакет программ на языке Delphi, проведены его отладка, тестирование и численные эксперименты.

3. РЕШЕНИЕ КРАЕВОЙ ЗАДАЧИ МЕТОДОМ БИХАРАКТЕРИСТИК

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

Здесь индекс $ * $ придается размерным величинам; индекс $m$ относится к материалу, в котором скорость продольных волн является наибольшей; $L{\text{*}}$ – характерный линейный размер;

$c_{1}^{{(k) * }} = \sqrt {\tfrac{{{{\lambda }^{{(k) * }}} + 2{{\mu }^{{(k) * }}}}}{{{{\rho }^{{(k){\text{*}}}}}}}} ,\quad c_{2}^{{(k) * }} = \sqrt {\tfrac{{{{\mu }^{{(k) * }}}}}{{{{\rho }^{{(k) * }}}}}} $

суть скорости распространения продольных и поперечных волн в $k$-й среде.

После введения безразмерных величин из уравнений (2.1), (2.2) после простых преобразований имеем ($i,j,k = 1,\;2$):

Уравнения (3.1) представляют собой линейную неоднородную гиперболическую систему дифференциальных уравнений первого порядка с постоянными коэффициентами. Ее характеристические поверхности в трехмерном пространстве (${{x}_{1}}$; ${{x}_{2}}$, $t$) представляют собой гиперконусы с осями, параллельными оси времени. Система уравнений (3.1) имеет два семейства характеристических конусов. Эти конусы совпадают с бихарактеристиками уравнений (3.1). Вдоль характеристик, лежащих в плоскости ${{x}_{\alpha }} = {\text{const}}$, уравнения (3.1) являются функциями только двух переменных (${{x}_{j}}$; $t$) ($j \ne \alpha $). Это обстоятельство указывает на то, что условия на бихарактеристиках могут быть получены как условия на характеристиках в соответствующей одномерной задаче. Соответствующие преобразования можно выполнить, если в системе уравнений (3.1) поочередно зафиксировать одну из пространственных переменных. При этом система уравнений (3.1) расщепляется на две системы уравнений, соответствующие направлениям $j = 1$ и $j = 2$ ($i = 1,\;2$): где введены обозначения Дифференциальные уравнения характеристик имеют вид ($i,j = 1,\;2$): а условиями на бихарактеристиках являются следующие: Здесь $\lambda _{{ij}}^{{(k)}} = c_{1}^{{(k)}}{{\delta }_{{ij}}} + c_{2}^{{(k)}}(1 - {{\delta }_{{ij}}}).$ Из (3.4) видно, что на каждой из двух гиперплоскостей имеются две пары семейств характеристик, определяющих продольные $c_{1}^{{(k)}}$ и сдвиговые $c_{2}^{{(k)}}$ скорости распространения волн. В каждой из двух плоскостей (${{x}_{j}}$, $t$) имеются по два семейства бихарактеристик положительного и отрицательного направлений. Верхний знак соответствует характеристикам положительного, а нижний – отрицательного направлений. Уравнения (3.4) и (3.5) соответствуют друг другу при одинаковой паре индексов и при одинаковом расположении знаков. Уравнения (3.2) и условия (3.5) используются для отыскания решения задачи [1]–[7].

Расчетная схема среды и штампа такова, что допускает систему координат ${{x}_{i}}$ ($i = 1,\;2$), в которой граничные поверхности являются координатными. Для построения алгоритма ${{D}_{1}}$ и ${{D}_{2}}$ разбиваются на ячейки, образуемые пересечениями координатных поверхностей ${{x}_{i}} = {\text{const}}$ ($i = 1,\;2$). Линейные размеры этих ячеек в направлении осей ${{x}_{1}}$ и ${{x}_{2}}$ считаются равномерными и равными $h$. Пересечения линий ${{x}_{i}} = {\text{const}}$ ($i = 1,\;2$) образуют узлы. В этих узловых точках отыскиваются значения искомых функций ${v}_{i}^{{\left( k \right)}}$, $\sigma _{{ij}}^{{\left( k \right)}}$ ($i,j = 1,\;2$) в различные моменты времени ${{t}_{n}} - \tau ,\;{{t}_{n}},\;{{t}_{n}} + \tau $ ($n = 1,2,\; \ldots ,\;N$) с шагом по времени $\tau $. Получившаяся сетка является трехмерной. Точечная сетка, на основе которых строится разностная схема, помимо упомянутых узловых точек, содержит точки, образованные пересечениями бихарактеристик с гиперплоскостями $t = {\text{const}}$. Использование явной разностной схемы второго порядка точности позволяет установить значения неизвестных величин в узловых точках плоскости ${{t}_{{n + 1}}}$ временного слоя по известным их значениям в узлах предыдущего слоя ${{t}_{n}}$. Принимается шаблон, состоящий из узла О и точек $E_{{ij}}^{{ \pm (k)}}$, лежащих на координатных линиях ${{x}_{j}} = {\text{const}}$ и отстоящих от точки О на расстояния $\lambda _{{ij}}^{{(k)}}\tau $. Наклонные прямые, исходящие из точки А, являются бихарактеристиками. В дальнейшем значениям функций в точке О приписывается верхний знак “О”; в точках $E_{{ij}}^{{ \pm (k)}}$ – нижний ($ij$) и верхний знак $ \pm $ (например, $\sigma _{{ij}}^{{ \pm \left( k \right)}}$), а в точке А дополнительный индекс не приписывается. Точки О и А представляют собой одну и ту же точку тела в моменты времени, отстоящие друг от друга на один шаг $\tau $ по времени.

На основании описанных точечных схем разрабатываемая ниже методика решения динамических задач позволяет определить скорости частиц ${v}_{i}^{{\left( k \right)}}$ и компоненты тензора напряжений $\sigma _{{ij}}^{{\left( k \right)}}$ в точке А на нечетном слое по времен ${{t}_{n}}$, если известны их значения на предыдущем слое ${{t}_{{n - 1}}}$ ($n = 1,\;2,\; \ldots ,\;N$) в точке О и в прилегающих к ней точках $E_{{ij}}^{{ \pm (k)}}$. Разностные схемы такого типа являются явными и удобны тем, что эти системы решаются последовательно от одного временного слоя к следующему. При этом искомые величины в каждой узловой точке, в отличие от неявной разностной схемы, вычисляются независимо от других.

4. МОДЕЛИРОВАНИЕ ДИНАМИКИ УПРУГОГО МАССИВА ПРИ ОБРАЗОВАНИИ ГОРИЗОНТАЛЬНЫХ ТРЕЩИН РАЗРЫВА

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

Здесь моделируется горизонтальная трещина разрыва, параллельная дневной поверхности. Расчеты произведены для среды с безразмерными параметрами: $\rho = 1$, ${{c}_{1}} = \sqrt {{{\gamma }_{{11}}}} = 1$, ${{c}_{2}} = \sqrt {{{\gamma }_{{12}}}} = 0.577$ при глубине и ширине трещины соответственно $L = 1$, $d = 0.1$, на интервале времени (0, 6).

Скачок напряжений на трещине задается в виде вертикального импульса (см. фиг. 2):

В расчетах шаги по пространственной сетке ${{h}_{1}} = {{h}_{2}} = h = 0.05$, по времени – $\tau = 0.025$, a параметр дельтаобразной функции $\varepsilon = h$.

5. ДИНАМИКА СРЕДЫ ПРИ СБРОСЕ ГОРИЗОНТАЛЬНЫХ НАПРЯЖЕНИЙ НА ТРЕЩИНЕ

При образовании сдвиговой трещины происходит сброс горизонтальных напряжений. На фиг. 10, 11 представлены осциллограммы скоростей и напряжений в эпицентре и фиксированных точках поверхности и среды. Скорости распространения продольных и поперечных волн в безразмерных величинах составляют ${{c}_{1}} = 1$ и ${{c}_{2}} = 0.577$ соответственно. Отсюда следует, что объемные (продольные) и сдвиговые (поперечные) волны от трещины к точке О(0, 0) придут за время ${{t}_{1}} = 1$ и ${{t}_{2}} = 73$ соответственно, а в точку А(0, 1) за ${{t}_{1}} = 35$ и ${{t}_{2}} = 2.33$, в точку В(0, 2) за ${{t}_{1}} = 2.15$ и ${{t}_{2}} = 3.7$ и в точку С(0, 3) за ${{t}_{1}} = 3.07$ и ${{t}_{2}} = 5.3$. В указанные моменты времени возникают колебания дневной поверхности. Характер колебаний в эпицентре на поверхности резко отличается от колебаний при образовании трещины отрыва. Преобладают горизонтальные колебания поверхности.

Скорость ${{{v}}_{1}} = 0$ в эпицентре О (и по всей прямой ${{x}_{2}} = 0$) объясняется типом воздействия на трещине. Так как в перпендикулярном направлении к дневной поверхности движется сдвиговая волна, то можно заметить, что чем ближе точка к эпицентру, тем выше амплитуда колебания горизонтальных скоростей точек дневной поверхности. Отрицательное направление скорости ${{{v}}_{2}}$ в точке О в промежуток времени (1; 7) объясняется взаимодействием фронтов распространяющихся волн.

На фиг. 12 представлены изолинии первого и второго инвариантов тензора напряжений в моменты времени $t = 1.5$ при образовании трещин сдвига. Они характеризуют распределение давления и интенсивности касательных напряжений в окрестности очага землетрясения.

На фиг. 13 представлено векторное поле скоростей в среде при сбросе горизонтальных напряжений на трещине:

Поскольку картина антисимметрична относительно вертикальной оси, то далее на фигурах представлена только правая половина (в левой половине направление векторов меняется на противоположное). В момент времени $t = 1$ объемная волна сжатия-расширения пробегает расстояние до поверхности, а сдвиговая волна 0.577$L$. Можно заметить вихревую зону, распространяющуюся вверх и вниз от трещины, порождаемую сдвиговыми деформациями на трещине. Ее окаймляет и опережает зона объемных расширений – сжатий, интенсивность скоростей перемещений в которой меньше, чем за фронтом сдвиговой волны.

Взаимодействие объемных и сдвиговых волн с дневной поверхностью и между собой (после отражения) порождают разные волны, среди них поверхностные. На фиг. 14 можно видеть, как начинает формироваться волна вдоль дневной поверхности, образуя поверхностную волну Рэлея. В данной среде она распространяется со скоростью ${{c}_{R}} = 0.918{{c}_{2}}$, т.е. близко к скорости распространения сдвиговой волны. Волны, идущие от трещины вниз, не встречая препятствий, рассеиваются с глубиной и с течением времени.

6. ДИНАМИЧЕСКОЕ НАПРЯЖЕННОЕ СОСТОЯНИЕ УПРУГОГО ШТАМПА ПРИ ПРЕЛОМЛЕНИИ СЕЙСМИЧЕСКИХ ВОЛН

Исследовалось напряженно-деформированное состояние штампа при преломлении сейсмических волн, при образовании трещин отрыва и сдвига. Проводился сопоставительный анализ в зависимости от расстояния штампа до эпицентра. Здесь на фиг. 15–19 представлены результаты расчетов для штампа, расположенного в эпицентре (фиг. 15а–19а) и на некотором расстоянии от него (фиг. 15б–19б), при образовании трещины разрыва в опорном массиве.

В качестве материала для штампа взят бетон с параметрами: $\rho $ = 2500 кг/м³, ${{c}_{1}}$ = 3593 м/с, ${{c}_{2}}$ = = 2200 м/с, при следующих размерах штампа: ${{d}_{1}} = 1$, ${{d}_{2}} = 0.5$.

Здесь представлены расчеты напряженно-деформированного состояния штампа при образовании в подстилающем пространстве трещины разрыва. На фиг. 15 представлены векторные поля скоростей точек тела ${{D}_{2}}$ в момент времени, когда преломленные волны дошли до середины штампа. При ${{d}_{3}} = 0$ (фиг. 15а) распространяется только продольная волна, и можно заметить эффект взаимодействия с боковой поверхностью. А при ${{d}_{3}} = 5$ (фиг. 15б) за продольной волной следует и поперечная волна, что соответствует типу воздействия. Здесь тоже заметен эффект взаимодействия, но сильнее с правой стороной. Это объясняется тем, что штамп стоит справа от эпицентра.

На фиг. 16 представлены векторные поля скоростей точек тела ${{D}_{2}}$ в момент времени, когда преломленная волна только добежала до верхнего торца штампа. На первой фигуре можно заметить, что за продольной волной начинается образование слабых поперечных волн, а на втором видно, что отраженная с правой боковой стороны волна порождает поперечную волну, которая достигла левой стороны.

На фиг. 17 представлены векторные поля скоростей точек тела ${{D}_{2}}$ в момент времени, когда преломленные волны отразились от верхнего торца. Здесь наблюдается сложная дифракционная картина. На фиг. 17а видно, что верхние угловые точки работают как источники продольных и поперечных волн.

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

Сравнивая их, видим, что тело в эпицентре подвергается обеим деформациям схожего характера, причем они симметричны оси ${{x}_{1}}$, и для тела, смещенного от эпицентра, характер объемных и сдвиговых изменений схож, но здесь волна деформации распространяется и от правой боковой стороны. А распределение давления и интенсивности касательных напряжений в первом теле сравнительно больше, чем во втором. Это связано с тем, что распространяющаяся волна возмущений теряет энергию по мере прохождения пути. Эти фигуры наглядно демонстрируют, если посмотреть на движение боковых поверхностей штампа, как деформируются здания при землетрясениях, причем сдвиговые деформации наиболее разрушительны.

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

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

ЗАКЛЮЧЕНИЕ

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

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

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