Введение. В 70-х годах прошлого века началось широкое использование метода конечных элементов (МКЭ) при анализе статической устойчивости склонов и откосов [1–3]. Как показала практика геотехнических расчетов МКЭ в отличие от упрощенных методик позволил учесть такие важные факторы как реальную геометрию и послойную структуру исследуемых объектов, а также наличие противооползневых сооружений и физическую нелинейность материала с заранее неизвестной областью пластичности [4]. В настоящее время многие исследователи стали применять МКЭ для анализа динамической устойчивости склонов реального заложения, а также прогнозирования рисков обвалов горных пород в подземных выработках и со склонов, расположенных вдоль автомобильных и железнодорожных магистралей. Также одним из актуальных направлений горной динамики является конечно-элементное моделирование последствий землетрясений с учетом кинематической нестабильности конгломераций [5]. В настоящее время существуют следующие методы динамических расчетов в области геотехники:

1. SRM (strength reduction method) — метод снижения прочности [6][7]. Разработан для расчета запаса прочности горного массива в физически нелинейной постановке. В SRM фактические прочностные параметры грунта итерационно уменьшаются путем деления на некоторый коэффициент больше 1:

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

2. LEM (limit equilibrium method) — метод предельной кинетики, базирующийся на принципе Даламбера [8][9]. Ориентирован на анализ динамической устойчивости сочлененных массивов горных пород.

3. TLEM (thin layer element method of FEM) — метод тонкослойных конечных элементов [10], в котором упругопластические элементы тонкого слоя используются для моделирования поведения кинематически нестабильных структур.

Анализ результатов, полученных с помощью методов SRM, LEM и TLEM, показал, что в настоящее время отсутствует единая концепция математического моделирования поведения структурно нестабильных геотехнических систем при нестационарном внешнем воздействии. Это обусловливает актуальность разработки методики динамического анализа систем типа «основание — ослабленный слой — блок» в конечно-элементной постановке с использованием нового подхода к моделированию плоскостей скольжения.

Материалы и методы. Уравнение движения механической системы в конечно-элементной формулировке представим в виде [11]:

(1)

где [M], [C], [K] — матрицы масс, демпфирования и жесткости ансамбля конечных элементов соответственно;

— векторы-столбцы соответственно узловых ускорений, скоростей, перемещений;  — векторы-столбцы заданных статических и динамических нагрузок соответственно в момент времени t . В дальнейшем полагаем, что матрицы [M] и [K] согласованные.

Для численного интегрирования уравнения (1) используем метод Ньюмарка [12]. При этом шаг интегрирования по временной оси t назначаем так, чтобы с достаточной точностью учитывались вклады физически значимых собственных пар. В дальнейшем рассмотрим кинематические способы возбуждения колебаний, заданные с помощью либо модельной сейсмограммы , либо с помощью модельной акселелограммы . При таком способе задания динамического воздействия второе слагаемое правой части уравнения (1) будет равно нулю:

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

Функцию представим в виде [13]:

, (2)

где A — начальная амплитуда;  — коэффициент затухания; — угловая частота внешнего воздействия. На рис. 1. приведен график функции для значений: A = 0,01553 м; = 0,7143; = 5 с–1.

Результаты исследования. В качестве первого модельного примера рассмотрим задачу о вынужденных колебаниях склона, расположенного на основании (рис. 2). Граничные условия задачи показаны на рис. 3, где буквой S обозначена исследуемая точка.

Соотношения геометрических параметров склона и основания (рис. 2) представлены в таблице 1.

Таблица 1

Соотношения геометрических параметров

Механические характеристики материала склона и основания: модуль деформации Е = 21 МПа; коэффициент Пуассона v = 0,3; удельный вес y = 1702 кг/м3; сцепление с = 45 кПа; угол внутреннего трения ф = 15°.

Для моделирования склона и основания используем объемные конечные элементы SOLID45 программного комплекса ANSYS Mechanical. Конечно-элементная модель для варианта с параметрами: ls = 2Hs, Hs = 10 м, отнесенная к глобальной декартовой системе координат, приведена на рис. 4.

Конечно-элементная сетка построена таким образом, что на поверхности контакта смежные узлы основания и склона имеют одинаковые координаты, но разные номера. Это принято для того, чтобы разместить в этом месте ослабленный слой. Кинематическое воздействие в форме модельной сейсмограммы (2) задаем на каждом шаге интегрирования ti в виде узловых перемещений на торцовых поверхностях модели с параметрами: X=0 и X=l1 + ls + l2.

Моделирование ослабленного слоя (рис. 4) выполняем с помощью упруго-вязких комбинированных конечных элементов COMBIN14 [14]. Двухузловой элемент COMBIN14, состоящий из пружины жесткостью k и демпфера жидкостного трения с коэффициентом демпфирования сv, приведен на рис. 5. В рассматриваемом случае этот элемент работает только на растяжение-сжатие.

В каждом узле контактной поверхности (рис. 4) вдоль глобальных осей X, Y, Z вводим элементы COMBIN14. Параметры комбинированных элементов:

k = 30 кН/м; ky = kz = 9,44∙10 7 кН/м; сv = 0,5.

В этом примере и далее вводим допущение о естественном недеформированном состоянии системы «основание — ослабленный слой — склон». Для вычислений используем нелинейный решатель комплекса ANSYS Mechanical.

Результаты конечно-элементного моделирования в виде визуализации деформированного состояния системы «основание — склон» с учетом максимального горизонтального смещения и распределения амплитудных горизонтальных перемещений Wx(t) показаны на рис. 6 и 7. Шаг интегрирования уравнения (1) t = 0,01 с. Как видно, введение 3D упруго-вязких элементов позволяет смоделировать эффект кинематической нестабильности механической системы «основание — ослабленный слой — склон» при кинематическом способе возбуждения колебаний.

Амплитудное значение перемещения в точке S составило WxS max  = 1,7 см. Для варианта склона ls = Hs (рис. 2) WxS max = 1,1 см. Графики колебаний основания и склона в исследуемой точке S (рис. 3) внаправлении оси X приведены на рис. 8.

На основании приведенных графиков видно, что начиная с момента времени t > 1,5 с наблюдается рассогласование колебаний основания и склона.

Рассмотрим поведение системы «основание — ослабленный слой — склон» (рис. 3) при возбуждении колебаний с помощью модельной акселерограммы . С этой целью дважды продифференцируем выражение (2). В результате получим:

(3)

График функции (3) для параметров: A = 0,01553 м; X = 0,7143; = 5 с–1 показан на рис. 9.

Кинематическое воздействие в форме модельной акселерограммы (3) по аналогии с сейсмограммой (2) задаем на каждом шаге интегрирования ti в виде узловых ускорений на торцовых поверхностях модели с параметрами: X=0 и X=l1 + ls + l2. На рис. 10 представлены графики колебаний в исследуемой точке S (рис. 3) при кинематическом воздействии в форме модельной акселерограммы. Сравнивая графики колебаний, приведенные на рис. 8 и 10, устанавливаем, что они практически совпадают. Это свидетельствует о корректности разработанной конечно-элементной модели, позволяющей описывать поведение системы «основание — ослабленный слой — склон» при различных способах нестационарного кинематического воздействия.

В качестве второго модельного примера рассмотрим задачу о вынужденных колебаниях склона с кинематически нестабильным клиновидным включением (рис. 11). В связи с симметрией конфигурации, в расчетной схеме учитываем только 1/2 часть склона и включения. Граничные условия для принятой расчетной схемы показаны на рис. 12. Здесь буквой S обозначена исследуемая точка, принадлежащая одновременно основанию склона и клиновидному включению.

Конечно-элементная модель склона и клиновидного включения показана на рис. 13. Как и в предыдущем примере, в данном случае используем элементы SOLID45 и COMBIN14 с теми же характеристиками материала.

Функция, описывающая модельную акселелограмму, имеет вид:

где A — амплитуда ускорения; — частота внешнего воздействия. На рис. 14 приведен график при A = 2,5 м/с2, = 2 Гц. Значения ускорения на i-том шаге интегрирования уравнения движения (1) прикладываем к узлам поверхности модели с координатой X = 0 (рис. 12).

Результат моделирования в виде распределения амплитудных значений перемещений Wx (t) показан на рис. 15. Шаг интегрирования t = 0,01 с. Графики колебаний основания склона и клиновидного включения в исследуемой точке S (рис. 12) в направлении оси X приведены на рис. 16.

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

«Дрейф»  на рис. 16 обусловлен тем, что данная конечно-элементная модель не имеет связей, препятствующих смещениям вдоль оси X. Как показано в [11], решить проблему «дрейфа» можно путем вычитания из значений перемещений  и  перемещения основания склона, которое представляет собой смещение «как жесткое целое». Отметим, что полученные амплитудные значения перемещений  и  позволяют оценить динамические параметры системы «склон — ослабленный слой — клиновидное включение».

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