← к сайту

26–30 мая 2026 года в Санкт-Петербурге прошёл XIII Международный научно-практический форум «Инновационные направления в проектировании горнодобывающих предприятий: безопасное и эффективное освоение месторождений полезных ископаемых». На форуме был представлен доклад о выполненной нами разработке и тестировании трёхмерной FDEM-модели для оценки устойчивости уступов и бортов карьеров.

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

Почему потребовалась трёхмерная FDEM-постановка

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

3D FDEM решает другую задачу. Метод конечных-дискретных элементов объединяет расчёт деформаций и напряжений с возможностью развития разрушения. В одной модели можно получить несколько механизмов разного масштаба: локальный вывал уступа, разрушение группы уступов и даже всего борта. Кроме того благодаря эксплицитной постановке возможно определение запаса устойчивости для всего карьера с учетом потенциальных обрушений любого масштаба.

Расчётная программа: FEM, DEM и ускорение на GPU

Разработанная программа использует идею FDEM: упругая работа породы описывается конечными элементами, а разрушение — интерфейсными элементами между ними. Объёмная часть модели представлена тетраэдрами, интерфейсы отвечают за раскрытие, сдвиг, разупрочнение и переход к остаточной прочности.

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

Расчёты выполняются на GPU с использованием CUDA. Целевой размер модели — более 1 млн элементов. Для диапазона 10⁴–10⁶ элементов GPU-расчёт может быть в десятки и сотни раз быстрее CPU. В программу включены модули импорта геомеханической модели, экскавации, принудительной стабилизации, порового давления и метода снижения прочности.

Плавная экскавация

Для карьерной модели важно корректно удалить элементы внутри проектной чаши. Если карьер мгновенно «вырезать» из массива, расчёт получит неестественную динамическую нагрузку, а часть разрушений может стать следствием искусственного "удара", а не реального напряжённого состояния.

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

Поровое давление и метод снижения прочности

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

Устойчивость оценивается методом снижения прочности SRM. Прочностные параметры массива уменьшаются пошагово коэффициентом SRF, начиная с 1,0. В отличие от ряда классических постановок снижается не только сдвиговая прочность, но и прочность на растяжение UTS. Коэффициент запаса устойчивости принимается равным критическому SRF в момент обрушения.

Преимущество трёхмерной FDEM-модели в том, что SRM можно выполнить по всей модели, а не только на заранее выбранных участках. Модель сама показывает, где возникнет критический механизм, когда будут превышены смещения и какие участки вовлекаются в разрушение. При этом возникновение разрушений на том или ином участке не мешает продолжить моделирование со снижением SRF далее, в отличие от имплицитного FEM. Благодаря этому в одной модели может быть выполнен расчёт устойчивости для любого множества участков.

Верификация отдельных алгоритмов

Перед расчётом реального карьера были проверены отдельные алгоритмы вычислительной программы-солвера. Поровое давление верифицировалось на простейшей схеме: между двумя конечными элементами задан интерфейсный когезионный (трещинный) элемент, в котором задается давление. При давлении Pв = 1,001 mg раскрытие между элементами должно происходить, при Pв = 0,999 mg — нет. Такой тест контролирует знак нормали, площадь грани и перевод давления в узловые силы.

Из приведенных выше видео можно увидеть, что алгоритм работает корректно с точностью до 0,1%.

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

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

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

Контактная модель проверялась на раскрытие, сдвиг, нелинейное допредельное поведение, запредельное разупрочнение и остаточную прочность. Для сдвига использовалось условие tau = C + sigma_n · tg(phi), а после потери сцепления — tau = sigma_n · tg(phi_r). Также проверялось удаление интерфейса при достижении предельного смещения.

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

Далее выполнено сравнение расчётов устойчивости простого однородного откоса. 3D FDEM дал коэффициент запаса устойчивости 1,59, расчёт в 2D Prorock — 1,67. Значения и форма разрушения сопоставимы, что подтверждает корректность перехода от тестов отдельных элементов к геомеханической задаче откоса.

Реальная 3D-модель карьера

Промышленный пример выполнен для проектного контура реального карьера на Дальнем Востоке РФ. Сетка сгущалась в областях, где ожидались потенциальные обрушения. Модель включала 329 211 объёмных элементов и 576 686 интерфейсных элементов. Средний размер элемента у контура составил около 5 м, минимальный показатель качества Scaled Jacobian — не ниже 0,20.

Геометрия строилась по STL-поверхностям: рельефу, разломам и границам литолого-структурной модели. Сетка была создана в Cubit. Габариты модели составили 1019 × 1013 × 327 м; высота бортов достигала 250 м, ширина карьера — до 725 м.

Свойства переносились из геомеханической блочной модели в конечные элементы по ближайшей ячейке, то есть по расстоянию между центроидами. В расчёт попадали код горной породы, GSI, UCS и производные параметры: сцепление, угол внутреннего трения, плотность, модуль деформации и другие показатели. Для зоны влияния буровзрывных работ использовалась глубина 50 м: параметр D равен 1,0 на поверхности карьера и убывает до 0 на глубине 50 м. Автоматизация переноса выполнена с применением PyVista, meshio и SciPy.

Структура массива, вода и напряжённое состояние

Системная трещиноватость задавалась по структурным доменам. Если грань интерфейсного элемента попадала в диапазон ориентации заданной системы трещин, ей назначались свойства трещины; в остальных случаях использовались свойства массива. Такой подход позволяет учесть не только прочностную неоднородность пород, но и пространственную ориентировку трещин.

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

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

Расчёт и последовательность разрушения

После экскавации при SRF = 1,0 все участки карьера оставались устойчивыми. Далее выполнялось снижение прочности, которое показало развитие разрушений нарастающего масштаба. Расчёт включал около 900 000 шагов и выполнялся на GPU RTX 4090. Базовый расчёт (стабилизация и экскавация) занял около 1 часа, этап снижения прочности — около 11 часов. Для трёхмерной модели такого размера это приемлемое время для практических задач, особенно с учётом того, что за раз анализируется весь карьер, а не один разрез.

В одном расчёте были выявлены четыре типа обрушений разного масштаба и генезиса. При SRF = 1,30 проявлялись клиновидные вывалы уступов на восточном борту, связанные с системной трещиноватостью. С повышением SRF сформировались более сложные механизмы, включая обрушения из-за трещин и отжима под весом восточного борта, уступы в слабых породах на юго-востоке, а также разрушение сложной формы большей части южного борта, ограниченное двумя крутонаклонными нарушениями и трещинами снизу.

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

Полная карта устойчивости вместо набора отдельных профилей

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

Карта показывает потенциальные обрушения любого типа и масштаба: уступов, групп уступов и борта. Такой результат отличается от классического FEM тем, что фиксируется не только участок с минимальным SRF, а коэффициент запаса определяется для любой точки на поверхности карьера.