CC BY
12
Вестник СПбГУ. Математика. Механика. Астрономия. 2023. Т. 10 (68). Вып. 2 УДК 519.63, 532.5 MSC 76U60, 76С05, 86A05
Численное моделирование
гидродинамических аварий:
размыв дамб и затопление территорий*
С. С. Храпов
Волгоградский государственный университет,
Российская Федерация, 400062, Волгоград, Университетский пр., 100
Для цитирования: Храпов С. С. Численное моделирование гидродинамических аварий: размыв дамб и затопление территорий // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2023. Т. 10(68). Вып. 2. С. 357-373. https://doi.org/10.21638/spbu01.2023.215
Построена математическая модель совместной динамики поверхностных вод и влекомых наносов, в которой учитываются нелинейная динамика жидкости и деформации дна. Динамика поверхностных вод описывается уравнениями Сен-Венана с учетом пространственно неоднородного распределения рельефа местности. Перенос влекомых наносов описывается уравнением Экснера, обобщенным на случай неоднородного распределения параметров подстилающей поверхности. Для численного интегрирования уравнений Сен-Венана и Экснера применяется устойчивый и хорошо апробированный метод CSPH-TVD второго порядка точности, параллельный CUDA-алгоритм которого реализован в виде программного комплекса «EcoGIS-Simulation» для высокопроизводительных вычислений на суперкомпьютерах с графическими сопроцессорами (GPU). Проведено гидродинамическое моделирование процессов размыва ограждающей дамбы реального гидротехнического объекта и затопления прилегающих территорий. Определены параметры прорана ограждающей дамбы и зон затопления, образовавшихся в результате развития гидротехнической аварии на хвостохранилище. На основе полученных результатов сделан вывод, что предложенный в работе метод численного моделирования совместной динамики поверхностных вод и влекомых наносов может являться более универсальным и эффективным (обладает существенно лучшей точностью и производительностью) по сравнению с существующими методиками расчета параметров прорана и зон затопления, особенно для гидродинамических течений со сложной геометрией на неоднородном рельефе дна.
Ключевые слова: математическое моделирование, уравнения Сен-Венана и Экснера, вычислительная гидродинамика, метод CSPH-TVD, зоны затопления, мелкая вода, влекомые наносы, прорыв и размыв дамб, хвостохранилище, параллельные вычисления, CUDA-алгоритм.
1. Введение. Деформация поверхности дна при слабых течениях сопровождается образованием структур в виде ряби, волн, рифелей, гряд, дюн и т.д., которые могут быть отслежены при исследовании донных отложений. В случае сильных течений происходит существенный размыв грунта, например при прорыве дамб, и образуется канал (проран), который быстро эволюционирует (углубляется и расши-
* Работа выполнена при финансовой поддержке Министерства науки и высшего образования РФ (госзадание №0633-2020-0003).
© Санкт-Петербургский государственный университет, 2023
ряется) под воздействием возрастающей интенсивности потока воды. К настоящему времени динамика поверхностных вод и влекомых наносов хорошо изучена как на основе экспериментальных данных, так и с использованием хорошо известных и апробированных математических моделей [1—8]. На основе этих моделей разработаны методики (см., например, [9, 10]) по расчету развития гидродинамических аварий на промышленных объектах, возникающих при прорыве ограждающих дамб. Данные методики используются как основа экспертной оценки безопасности при проектировании гидротехнических сооружений и определении границ зон затопления в окрестности промышленных объектов. Как правило, аналитические формулы, используемые в этих методиках, основаны на решениях упрощенных (линеаризованных и одномерных) уравнений, являющихся некоторыми приближениями исходной нелинейной двумерной модели динамики поверхностных вод и наносов, что, в свою очередь, приводит к увеличению погрешности расчетов, ограничению класса решаемых задач по геометрии водных объектов и структуре течения, а также требует больших временных затрат при выборе параметров моделей и нормировочных коэффициентов.
В настоящее время для изучения гидродинамических течений в природных и технических системах наряду с физическим экспериментом и аналитическими методами исследования математических моделей применяется важный и мощный инструмент — вычислительный эксперимент, основанный на численном решении уравнений гидродинамики хорошо апробированными методами и применении параллельных технологий для повышения производительности расчетов [11—18].
Цель данной работы — разработка эффективной методики расчета зон затопления и размыва грунта дамб гидротехнических объектов (хвостохранилищ) с произвольной геометрией на сложном рельефе местности, основанной на прямом численном интегрировании уравнений нелинейной динамики поверхностных вод и влекомых наносов. В разделе 2 описана математическая модель совместной динамики поверхностных вод и влекомых наносов, а в разделе 3 — численный алгоритм и структура параллельного вычислительного модуля программного комплекса «Еес^К-81ти1а1юп». Результаты численного моделирования гидродинамической аварии на реальном гидротехническом сооружении представлены в разделе 4.
2. Математическая модель динамики поверхностных вод и наносов.
Динамика поверхностных вод с учетом основных физических и метеорологических факторов, влияющих на процессы затопления территорий, описывается системой уравнений Сен-Венана:
где Ь — функция рельефа местности; Н — толщина слоя поверхностной воды; П = Н + Ь — уровень свободной водной поверхности; и = {и, V} — вектор скорости течения поверхностного слоя жидкости; = {д/дх, д/ду} — дифференциальный оператор набла в плоскости (х,у); д = 9.81 м/с2 — ускорение свободного падения. Величина ц определяет скорость притока/оттока жидкости за счет действия источников/стоков в поверхностном слое воды [19]. Уравнения Сен-Венана основаны на модели мелкой воды [21], в которой уравнения гидродинамики усредняются по вертикальной ^-координате в предположении однородной несжимаемой
(1)
жидкости, находящейся в гидростатическом равновесии. Условие гидростатического равновесия слоя жидкости означает, что мы пренебрегаем вертикальной компонентой скорости жидкости, считая ее среднее значение в слое равным нулю, а на свободной поверхности жидкости задается условие равенства давлений: p(r¡) = pa, где p — гидродинамическое давление в жидкости, pa — атмосферное давление. С учетом этих предположений из z-компоненты уравнения движения жидкости получаем: p(z) = pa + gg(n — z), где q = const — плотность жидкости.
В общем случае суммарная удельная сила f, действующая на слой жидкости в модели мелкой воды, включает силы придонного (f>) и вязкого турбулентного трения f), Кориолиса (fo), ветра (fw) (см., например, [11, 20]).
В рамках стандартной модели квадратичного трения жидкости о дно водоема величина f> представима в виде
fA = ~Au|u|, (3)
где Л = 2gn2M/Н4/3 — коэффициент гидравлического трения, пм — коэффициент Маннинга, определяющий уровень шероховатости дна.
Величина удельной силы вязкого турбулентного трения имеет вид
f„ = ¿V±T, (4)
где T — тензор вязких напряжений с компонентами:
Тхх = 2,Нд£, Тт=2„н\Тху = Тух = ,Н^ + ^-у (5)
Для коэффициента турбулентной вязкости v в (5) имеем
где С ~ 0.4 — эмпирическая постоянная; Ь2 — площадь расчетной ячейки. Удельная сила Кориолиса:
^ = (1п х,1п у) =2(-у Псо8©,иПз1п©), (7)
где П — угловая скорость вращения Земли; © — географическая широта. Удельная сила ветра, действующая на слой жидкости:
¡ги = - и) | ЛУ - и | , (8)
Ид
где да — плотность атмосферы; ' = (тх,ту) — скорость ветра; Са — параметр, зависящий от состояния водной поверхности (Са ~ 10~3).
Динамика наносов описывается уравнением Экснера [22]: (1-ф)дЬ/д1+У Л = 0, где величина расхода наносов Л в общем случае определяется нелинейной зависимостью от скорости потока поверхностных вод. На неоднородном дне величина Л должна также зависеть и от наклона дна Ь = 0 (см., например, работы [5, 7, 8,
23, 24]) — J = Jb — cj|Jb|V^b, где Jb — расход наносов на плоском дне, а cj « 1.55 — эмпирическая постоянная [24], зависящая от типа и состояния грунта. Таким образом, для моделирования динамики размыва грунта и перемещения влекомых наносов в донных отложениях получим следующее уравнение:
db
(1 - ф) — + V±Jb = cjV± {|Jb|V±6} , (9)
где ф — пористость донного грунта. Расход наносов на плоском дне будем определять по формуле Грасса [23]:
(aj u |u|m, |u| >vCr
Jb = in II/ , (10)
[0, |u| < Vcr
где aj и m — коэффициенты, зависящие от характеристик переносимого материала; vcr = cv ¿50 H1/6 — критическая (неразмывающая) скорость, определяемая по формуле Шамова [2], cv « 3.7-6 — эмпирическая постоянная [2], зависящая от типа и состояния грунта.
Стоящее в правой части (9) слагаемое диффузионного типа связано с дополнительным потоком материала грунта на наклонном дне (V±b = 0), которое обеспечивает, в частности, наблюдаемое обрушение стенок канала за счет поперечного к направлению скорости жидкости потока материала грунта. Во многих моделях при замыкании уравнения Экснера аналитическими соотношениями для расхода наносов параметр m = 2 (см., например, [5, 6, 8]). Величину Aj можно определить следующим образом [4, 25]:
, (11)
(s — 1)vgH ¿50
где s = ps/p — удельная плотность наносов; ¿50 — медианная крупность наносов; пм — коэффициент Маннинга [20].
В качестве начальных условий для системы уравнений (1)-(9) задаются: функция рельефа дна b(x, y,t = 0) = b0(x, y), построение которой рассмотрено в разделе 4; распределение воды H(x, y,t = 0) > 0 с постоянным уровнем свободной поверхности n(x,y,t = 0) = no(x,y) = const и u(x,y,t = 0) = 0 внутри хвостохранилищ, ограниченных дамбами, а за их пределами — H(x,y,t = 0) = u(x,y,t = 0) = 0. Кроме того, при моделировании гидротехнической аварии в начальный момент времени на ограждающих хвостохранилище дамбах вносится возмущение на функцию рельефа дна (b(x, y,t = 0) = b0(x, y) + b(x, y)) для создания начального относительно небольшого прорана (канала), в котором формируется интенсивный поток воды, разрушающий дамбу и приводящий к затоплению прилегающих территорий.
Гидротехнические объекты размещаются внутри прямоугольной расчетной области x € [xmin,xmax] и y € [ymin,ymax], на границах которой задаются условия свободного протекания: V±f (x, y,t)|r =0 , где f = {b, H, u}.
3. Численный метод и программная реализация. Для численного интегрирования системы уравнений (1)-(9) воспользуемся устойчивым и хорошо апробированным методом CSPH-TVD второго порядка точности, подробно описанным в работах автора [11, 20] при решении практических задач динамики поверхностных вод на реальном рельефе местности. С учетом транспорта донных наносов,
описываемых уравнением (9), необходимо обобщить численный алгоритм метода СSPH-TVD, дополнив вектор консервативных переменных и, потоки Ё и О, а также вектор источников и сил Ф из [11, 20] новыми величинами. Здесь ограничимся кратким описанием основных этапов СSPH-TVD-алгоритма на декартовой сетке ), уделив внимание только ключевым моментам, связанным с модификацией метода из-за учета совместной динамики поверхностных вод и влекомых наносов. Лагранжев этап СSPH-TVD-метода запишем в виде
Ui,j(t + г) = Uitj(t) + гФ^- (Uij(t + 0.5r))
(12)
где т — шаг интегрирования по времени, а для нахождения величины Ф в момент времени t + 0.5т применяется процедура предиктор-корректор (leapfrog):
U
/я\
H u
\b J
(
Ф
q
gHV±(H + b) + H f 1 CjV±{\Jb\V±b}
\
\1 - ф
Аппроксимация пространственных производных в (12) осуществляется в соответствии с модифицированным SPH-подходом [11].
Эйлеров этап СSPH-TVD-метода можно представить в следующем виде:
и¿¿(г + т) = и + - - - ^1-1)2,3 + ~ , (13)
где дробные пространственные индексы г ± 1/2 и ] ± 1/2 соответствуют границам ячеек расчетной сетки; Н = Дж = Ду — пространственный размер ячеек; Ё4±1/2^ и Ог¿±1/2 — средние на временном промежутке + г] значения потоков массы и импульса вдоль соответствующих направлений (ж, у),
/ Hu \ Hv
Hu2 Huv
F = Huv , G = Hv2
Jbx Jby
\1 - ф/
\1 - ф/
Вычисления потоков Ё±1/2 ^ и СО ¿±1/2 на эйлеровом этапе CSPH-TVD-метода осуществляются на основе модифицированного TVD-подхода [11, 26] и приближенных методов решения задачи Римана для уравнений мелкой воды (Лакса — Фридрихса — LF, Хартена —Лакса —Ван Лира — HLL, HLLC) [27, 28].
Условие устойчивости численного алгоритма совместной динамики поверхностных вод и наносов (12), (13), основанного на CSPH-TVD-методе, имеет вид
Km
h
2 Up
h h2 v~'2D
(14)
где 0 < К < 1 — число Куранта; vv = max|u|; va = max (|u| + \/gH); D = max |Jb| /(1 — ф), а функции max и min вычисляются в момент времени t для всей пространственной области моделирования. Третий аргумент функции min в (14)
обеспечивает устойчивость явной численной схемы для дифференциального уравнения параболического типа (9) и дополняет условие устойчивости метода CSPH-TVD [11] с учетом динамики наносов.
Отметим, что численный алгоритм совместной динамики поверхностных вод и наносов имеет тот же второй порядок точности, как и CSPH-TVD-метод для уравнений Сен-Венана (см. [11]), поскольку для аппроксимации пространственных производных в уравнении Экснера (9) применяются SPH- и TVD-подходы, обладающие точностью 0(h2), а использование для интегрирования по времени алгоритма предиктор-корректор (leapfrog) обеспечивает точность 0(т2).
При разработке расчетного модуля с рабочим названием «EcoGIS-Simulation», предназначенного для моделирования самосогласованной динамики поверхностных вод и наносов («CUDA-SWD» + «CUDA-SD»), использовалась технология параллельных вычислений CUDA для графических процессоров (GPU). Реализация расчетного модуля «EcoGIS-Simulation» на основе параллельной технологии CUDA позволяет эффективно использовать компьютерные платформы (CPU+GPU) и повысить производительность вычислений в сотни раз по сравнению с последовательной версией на CPU.
GPU
Рис. 1. Потоковая диаграмма для расчетного модуля «EcoGIS-Simulation». Показаны последовательность выполнения блоков CUDA-ядер на GPU и потоки данных между CPU и GPU.
На рис. 1 показана структура вычислительного модуля «EcoGIS-Simulation» в виде потоковой диаграммы, демонстрирующей последовательность выполнения следующих CUDA-модулей: SWD — вычисление временного шага (14) и расчет динамики поверхностных вод (1), (2) на основе CSPH-TVD-метода, подробное описание численного алгоритма представлено в работах [14, 17]; SD — расчет динамики наносов (9) на основе обобщенного численного алгоритма CSPH-TVD-метода (12)-(14).
4. Результаты моделирования гидродинамической аварии. Рассмотрим особенности применения построенной математической модели и разработанного программного комплекса для моделирования аварий на гидротехнических сооружениях, связанных с размывом дамб и затоплением прилегающих территорий. Для демонстрации имитационного моделирования гидродинамических аварий с использованием нашего программного комплекса было выбрано хвостохранилище Орловской обогатительной фабрики (ООФ), расположенной на равнине в отсутствии каких-либо русловых стоков. При решении задачи по прогнозированию динамики развития гидродинамических аварий необходимо: построить матрицы высот моделируемого участка местности на основе данных дистанционного зондирования Земли (ДЗЗ), топографических карт и чертежей территории хвостохранилища; провести моделирование гидродинамической аварии в результате частичного разрушения ограждающей дамбы хвостохранилища; определить размеры прорана и скорость потока в
нем; спрогнозировать развитие потенциальной гидродинамической аварии с определением зоны затопления и указанием их границ; построить карты распределения глубины воды и скорости течения в зонах затопления; провести оценку общей площади и протяженности затопления, оценить площади затопления населенных пунктов, промышленной и дорожной инфраструктуры.
Описание объекта исследования и исходные данные. Объект исследования — хвостохранилище Орловской обогатительной фабрики, которое расположено в 140 км северо-западнее от г. Усть-Каменогорск и на 4.5 км западнее промышленной площадки ООФ (рис. 2, о,).
Высотные отметки различных резервуаров (прудков) хвостохранилища ООФ: А) уровни верхнего яруса ограждающей дамбы 274.5-275 м, дна прудка — 273 м, у подножия хвостохранилища — 264 м, высота ограждающей дамбы относительно дна прудка 1.5-2 м; В) уровни верхнего яруса ограждающей дамбы 274.5-277 м, дна прудка — 273 м, у подножия хвостохранилища — 265 м, высота ограждающей дамбы относительно дна прудка 1.5-4 м; С) уровни верхнего яруса ограждающей дамбы 270 м, дна прудка — 266 м, у подножия хвостохранилища — 263 м, высота ограждающей дамбы относительно дна прудка 4 м; D) уровни верхнего яруса ограждающей дамбы 269.5-270 м, дна прудка — 264 м, у подножия хвостохранилища — 263 м, высота ограждающей дамбы относительно дна прудка 5.5-6 м.
Площадь и объем различных прудков хвостохранилища ООФ: А) на уровне 274 м — 848 тыс. м2 и 738 тыс. м3; В) на уровне 274 м — 281 тыс. м2 и 233 тыс. м3; С) на уровне 269.5 м — 201 тыс. м2 и 696 тыс. м3; D) на уровне 269 м — 201 тыс. м2 и 967 тыс. м3.
Построение цифровой модели местности. Построение цифровой модели рельефа для целей гидродинамического моделирования проводилось в несколько этапов на основе топографических планов и данных дистанционного зондирования Земли ^ИТМ). Сначала строился локальный рельеф хвостохранилища на основе про-
Рис. 2. Положение хвостохранилища ООФ на карте (а). Трехмерная модель (матрица высот) моделируемого участка местности (б). Латинскими буквами (А, В, С, Ц) отмечены различные резервуары хвостохранилища ООФ.
ектных планов, который в дальнейшем сшивался с глобальным рельефом (SRTM). Затем проводилась дополнительная обработка матриц высот для устранения зубчатости рельефа (применялись различные методы сглаживания) и адаптации их к использованию в расчетном модуле программного комплекса «EcoGIS-Simulation» (проводилась нарезка матриц и конвертация в обменный формат .grd). В результате дополнительной обработки была получена матрица высот 1200 х 1200 ячеек с размером элемента 5 м. Моделируемый участок охватывает территорию 6 х 6 км (Xmin = Ушт = 0 м, Хтах = ymax = 6000 м).
Таким образом, на основе методов геоинформационного моделирования построена итоговая цифровая модель рельефа (ЦМР) местности с учетом проектных решений по реконструкции хвостохранилища ООФ, представленная на рис. 2, б в виде трехмерной карты высот.
Сценарий гидродинамической аварии. 1. Начальное состояние и сценарий гидродинамической аварии. Начальный уровень воды в хвостохранилище ООФ перед моделированием составляет: для резервуара A — по = 274 м (максимальная глубина прудка 1 м); для резервуара B — п0 = 274 м (максимальная глубина прудка 1 м), для резервуара C — п0 = 269.5 м (максимальная глубина прудка 3.5 м); для резервуара D — по = 269 м (максимальная глубина прудка 5 м). Для задания величины b считаем, что в начальный момент времени в результате механического локального повреждения (повреждение экскаватором при проведении земляных работ или обвала дамбы) формируется канал (траншея) шириной ~ 5 м (поперек дамбы), длиной 10-20 м (от дна прудка до внешнего края дамбы) и глубиной 1—1.5 м от верхней отметки дамбы, формируя таким образом канал с водным потоком прорыва глубиной порядка метра. Далее в результате воздействия потока воды происходят размыв дна канала и обрушение боковых стенок канала с образованием естественного откоса. После образования прорана селевой поток (смесь воды и грунта) устремляется к подножию дамбы и впоследствии образует зону затопления в окрестности хвосто-хранилища ООФ. При моделировании динамики затопления рассмотрены несколько вариантов (моделей) формирования прорана на различных участках дамб резервуаров A, B, C, D хвостохранилища ООФ: модель A — на восточной части дамбы резервуара A; модель B — на южной части дамбы резервуара B; модель C —на южной части дамбы резервуара C; модель D — на западной части дамбы резервуара D.
2. Параметры дамбы для моделирования размыва и затопления. Верхние ярусы дамб, которые будут размываться, состоят из насыпных грунтов (хвостов). Параметры грунтов, используемые в моделировании размыва дамбы и затопления, следующие: медианная крупность ¿50, пористость ф, плотность частиц грунта р. По данным инженерно-геологических исследований для основной фракции дамбы имеем: ф = 0.5, р = 3.8 г/см3, ä50 = 0.031 мм.
Средняя объемная доля первой фракции на верхних ярусах дамбы порядка 75 %. В результате механического воздействия и последующего обрушения грунта по бокам канала основная фракция может перемешиваться с другими более крупными фракциями (песок, щебень), среднее значение медианной крупности которых ¿50 — 1 мм. В численных расчетах для полученной смеси фракций при нормальных условиях примем ¿50 = 0.273 мм, пм = 0.02, cj = 1.5 и cv = 4.4 (см., например, [24, 25]).
Расчет образования прорана. Численное моделирование процесса разрушения дамбы в результате развития гидродинамической аварии на хвостохранилище ООФ
Модель А Модель В
Рис. 3. Структура прорана дамбы хвостохранилища ООФ на конечной стадии развития гидродинамической аварии в различных моделях. Положение контрольных точек показано красным символом.
проводилось на основе алгоритма самосогласованной динамики затопления и размыва грунта (12)—(14). Время моделирования составляло 24 ч от момента возникновения аварии. За это время резервуары A, B, C, D хвостохранилища ООФ полностью опустошаются. Результаты расчетов для моделей A, B, C, D приведены на рис. 3-4.
Структура сформировавшихся в процессе развития гидродинамической аварии на хвостохранилище ООФ проранов показана на рис. 3. Из рисунка видно, что поток воды и размывает дно канала, созданного в результате механического локального повреждения, и приводит к обрушению боковых стенок канала, увеличивая также и его ширину. На рис. 3 хорошо заметны намывные структуры внутри хвостохрани-лища и снаружи, которые возникают в результате оседания и переноса грунта дамб. Ширина образовавшегося в процессе размыва дамб прорана составляет 25-50 м: A — 30-35 м, B — 25-30 м, C и D — 45-50 м.
С ростом размера прорана увеличиваются скорость и глубина образовавшегося селевого потока (смеси воды и грунта), что в свою очередь приводит к еще большему размыву грунта в области прорана и быстрому увеличению объемного расхода воды. Временная зависимость объемного расхода воды Q(t), вытекающего через проран ограждающей дамбы хвостохранилища, показана на рис. 4, a. Для моделей A и B величина Q достигает максимального значения ~ 40-50 м3/с через 100-120 мин после прорыва дамбы. Кроме того, в модели A на временах 320-340 мин хорошо за-
500-1
""1""1""1""1""1Т;"1||"1||||1||||1
200 400 600 800 1000 t, мин
276-
272-
2 268-
N
264-
260-
1/"
11 11 11 11 11 I I 11 11 11 11 11 11 11 11 11 11 11 11 I И | 11 11111 П |
0 200 400 600 800 1000 t, мин
200
■ I■■II!■ II || ■ 1111
400 600 800 1000 t, мин
.I.■.■!■... I.... I
400 600 800 1000 t, мин
Рис.4- Характеристики прорана дамбы в зависимости от времени для различных моделей. Показаны: о — объемный расход воды, вытекающей через проран дамбы; уровень дна б, скорость в и глубина потока г в контрольной точке (см. рис. 3).
метен второй максимум со значением Qmax ~ 25 м3/с, который связан с размывом канала внутренней дамбы, разделяющей резервуар А на две части, и увеличением сброса за счет быстрого притока воды из внутреннего резервуара. Модели С и В характеризуются более быстрым размывом дамб и существенно большей величиной объемного расхода воды через образовавшийся проран, что обусловлено большей глубиной воды в прудке и геометрической структурой дамб, которые имеют меньшую ширину и больший уклон. Максимальный уровень сброса (расхода) воды со значениями Qmax ~ 350-500 м3/с достигается к моменту времени Ь « 40 мин в моделях С и В соответственно.
На рис. 4, б-г показаны временные зависимости параметров потока в окрестности контрольных точек, отмеченных на рис. 3 красными крестиками. В моделях А и В происходит плавный процесс размыва грунта и за время Ьр « 120-140 мин уровень дна канала в проране опускается на 4 м до отметки 269 м (рис. 4, б). На временах Ь > Ьр вода продолжает вытекать из прудков не размывая грунт, так как скорость потока оказывается меньше критической (неразмывающей) скорости усг (рис. 4, в, г).
Максимальные значения скорости и глубины потока для моделей Л и Б составляют 5 м/с и 0.5-1 м соответственно. В моделях С и Б происходит более интенсивный размыв грунта с быстрым уменьшением уровня дна канала в проране за ~ 40 мин. Для модели С уровень дна канала в проране опускается на 6 м, а в Б — на 7.5 м. На временах 40 мин < Ь < « 100 мин уровень дна канала в проране в этих моделях увеличивается на 1.5 м за счет обрушения боковых стенок канала. Максимальные значения скорости и глубины потока для моделей С и Б составляют 7 м/с и 3.3-4 м соответственно.
Расчет зон затопления. В процессе размыва ограждающей дамбы хвостохра-нилища ООФ формируется широкий и глубокий канал (проран), через который образовавшийся селевой поток устремляется к подножию дамбы и приводит к затоплению прилегающих территорий. Результаты гидродинамического моделирования и расчета зон затопления для различных моделей представлены на рис. 5-6.
На рис. 5 показаны временные зависимости объема воды в хвостохранилище ООФ V(Ь) и площади затопления прилегающих к дамбам территорий. В моделях Б, С, Б полное опустошение прудков хвостохранилища ООФ происходит менее чем за 4 ч. Для модели Л через 4 ч после начала прорыва объем воды в хвостохра-нилище уменьшается приблизительно наполовину, а через 16 ч в хвостохранилище остается менее 5 тыс. м3, при этом глубина воды в хвостохранилище оказывается меньше 5 см. Как видно из рис. 5, б, площади затопления в моделях С, Б практически совпадают и составляют ~ 4 км2. Несмотря на то что объем воды в прудке Л почти в три раза больше, чем в прудке Б, площадь затопления в модели Л оказывается приблизительно в 2.5 раза меньше, чем в модели Б. Это обусловлено особенностями рельефа прилегающих к дамбам территорий. На восточной стороне имеется ограниченная с севера дорожной насыпью котловина, в которую и выливается весь объем воды в модели Л, а с южной и западной сторон хвостохранилища рельеф приблизительно равнинный и не ограничен насыпями.
1111 Г] 11 ■ ■ I ■ ■ ч 111111 и 1111111111111111111111111111
0 200 400 600 800 1000
5x10 -а
4x10 Н
3x10°-
2x10 Н
10 0
Л /;
I.
1;>
гтттттттттттттттттттттттттттттттттттттттттттл
0 200 400 600 800 1000 t. мин
Рис. 5. Зависимость объема воды от времени в прудке хвостохранилища (а) и площадь затопления (б) для различных моделей.
t. мин
В технических приложениях границы зон потенциального затопления территорий определяются по максимальному значению глубины воды в расчетных ячейках за все время моделирования. Расчетные границы зон затопления прилегающих к
хвостохранилищу ООФ территорий для различных моделей (А, В, С, В) показаны на рис. 6. Распределение глубин в зонах затопления отображается различным полупрозрачным цветом. Минимальное значение глубины воды, по которому строилась граница затопления, составляет 0.05 м (цвет границы красный).
Общая характеристика зон затопления для различных расчетных моделей:
Модель А — затапливается только территория, прилегающая к хвостохранили-щу с восточной стороны. Зона затопления ограничена на севере дорожной насыпью. Максимальная площадь затопления составляет 1.2 км2, максимальная глубина воды — 2.5 м, максимальная скорость течения — 2.8 м/с.
Модель В — затапливаются территории, прилегающие к ООФ с южной и западной сторон. Максимальная площадь затопления составляет 2.5 км2, максимальная глубина воды — 1.3 м, максимальная скорость течения — 1.7 м/с.
Модель С — также затапливаются территории, прилегающие к ООФ с южной и западной сторон. Максимальная площадь затопления составляет 4.1 км2, максимальная глубина воды — 2.3 м, максимальная скорость течения — 5.7 м/с.
Модель В — затапливается только территория, прилегающая к хвостохрани-лищу с западной стороны. Максимальная площадь затопления составляет 4.2 км2, максимальная глубина воды — 2.2 м, максимальная скорость течения — 7 м/с.
В зонах затопления для моделей А, В, С, В отсутствует промышленная и жилая инфраструктура. Затапливаются только дороги протяженностью 1.5, 1.4 и 0.8 км в моделях А, В, С соответственно.
5. Выводы. В заключение сформулируем основные результаты работы и обсудим области их возможного практического применения.
Построенная математическая модель самосогласованной динамики поверхностных вод и влекомых наносов позволяет рассматривать гидродинамические течения в окрестности водоемов произвольной геометрии с учетом неоднородного распределения параметров подстилающей поверхности (донных отложений и грунтов), различные типы течений в селевых потоках (смеси воды и грунта) при размыве грунта ограждающих дамб, включая как продольный перенос наносов, так и поперечный, связанный с обрушением боковых стенок канала в образовавшемся проране. Модель основана на уравнениях Сен-Венана, описывающих динамику поверхностных вод в приближении теории мелкой воды, и уравнении Экснера, описывающего динамику влекомых наносов на наклонном дне. Представленный в работе численный алгоритм (лагранжево-эйлеров метод CSPH-TVD) решения задачи самосогласованной динамики поверхностных вод и влекомых наносов является условно устойчивым для широкого диапазона возможных скоростей потока, имеет второй порядок точности 0(т2 ,h2) и является двухшаговым типа предиктор-корректор. Устойчивость алгоритма, т.е. отсутствие осцилляций, в численном решении достигается благодаря применению TVD-подхода и приближенных методов решения задачи Римана при вычислении гидродинамических потоков. Существуют и другие устойчивые численные методы решения уравнений типа Экснера, имеющие второй порядок точности [16, 29, 30]. На основе построенной численной модели разработан параллельный CUDA-алгоритм, реализованный в виде программного комплекса для суперкомпьютеров с GPU с рабочим названием «EcoGIS-Simulation», и проведено гидродинамическое моделирование размыва грунта дамб и затопления прилегающих территорий при возникновении аварийных ситуаций на реальном гидротехническом объекте — хвостохранилище ООФ. Определены размеры, глубина и скорость потока в образовавшемся проране ограждающей дамбы. Построены границы зон затопления прилегающих к хвостохранилищу территорий, определены общая площадь и протяженность затопления промышленной и дорожной инфраструктуры.
Отметим, что методики расчета образования прорана в процессе размыва и разрушения дамб (см., например, [9, 10]), основанные на аналитических решениях упрощенных (линеаризованных и одномерных) уравнений динамики поверхностных вод и влекомых наносов, в общем случае, т. е. для реальной местности со сложной геометрией водных объектов и неоднородного распределения параметров подстилающей поверхности, являются существенно менее точными и более затратными по сравнению с предложенной в работе универсальной методикой численного моделирования самосогласованной динамики поверхностных вод и влекомых наносов. С учетом неоднородного распределения параметров реальных грунтов, сложной геометрии водных объектов, нестационарности потока воды и нелинейных эффектов, возникающих при взаимодействии гидродинамических потоков поверхностных вод и влекомых наносов, погрешность упрощенного подхода [10] может существенно возрасти и привести к некорректным (неадекватным) оценкам параметров прорана и границ зон затопления. Кроме того, применение в наших численных моделях высокопроизводительных параллельных вычислений на суперкомпьютерах с GPU позволяет эффективно рассчитывать динамику зон затопления и размыва подстилающей поверхности (грунта) на обширных территориях с высоким пространственным разрешением.
Для определения точности и пределов применимости рассматриваемой в работе модели самосогласованной динамики поверхностных вод и влекомых наносов требуется проведение как натурных экспериментов с реальными грунтами, так и модельных численных расчетов на основе полной трехмерной гидродинамики многослойных сред.
Благодарность. Автор признателен А. В. Хоперскову за обсуждение результатов работы и полезные замечания.
Литература
1. Леви И. И. Динамика русловых процессов. Ленинград, Госэнергоиздат (1957).
2. Караушев А. В. Речная гидравлика. Ленинград, Госэнергоиздат (1969).
3. Richards K. J., Taylor P. A. A numerical model of flow over sand waves in water of finite depth. Geophys. J. Int. 65, 103-128 (1981).
4. Van Rijn L. C. Sediment transport, Part I: Bed load transport. Journal of Hydraulic Engineering 110 (10), 1431-1456 (1984).
5. Petrov P. G. Motion of a bed load. J. Appl. Mech. Tech. Phys. 32 (5), 717-721 (1991).
6. Venditti J. G., Church M. A., Bennett S. J. Bed form initiation from a flat sand bed. J. Geophys. Res.: Earth Surface 110 (F01009), 1-19 (2005).
7. Alan G.D., Peter D.Th. Advances in the study of moving sediments and evolving seabeds. Surveys in Geophysics January 29 (1), 1-36 (2008).
8. Petrov A. G., Potapov I. I. Sediment transport under normal and tangential bottom stresses with the bottom slope taken into account. J. Appl. Mech. Tech. Phys. 55 (5), 812-817 (2014).
9. Временные методические рекомендации по расчету зон при внезапном прорыве ограждающих дамб хвостохранилищ. Белгород, ВИОГЕМ (1981).
10. Безопасность гидротехнических сооружений на объектах промышленности и энергетики: сборник документов. Москва, ЗАО НТЦ ПБ (2010).
11. Храпов С. С., Хоперсков А. В., Кузьмин Н.М., Писарев А. В., Кобелев И. А. Численная схема для моделирования динамики поверхностных вод на основе комбинированного SPH-TVD-подхода. Вычислительные методы и программирование: новые вычислительные технологии 12 (1), 282-297 (2011).
12. Bulatov O. V., Elizarova T.G. Regularized shallow water equations and an efficient method for numerical simulation of shallow water flows. Comput. Math. and Math. Phys. 51, 160-173 (2011).
13. Сухинов А. И., Чистяков А. Е., Семенякина А. А., Никитина А. В. Параллельная реализация задач транспорта веществ и восстановления донной поверхности на основе схем повышенного порядка точности. Вычислительные методы и программирование: новые вычислительные технологии 16 (2), 256-267 (2015).
14. Dyakonova T., Khoperskov A., Khrapov S. Numerical model of shallow water: the use of NVIDIA CUDA Graphics processors. Communications in Computer and Information Science 687, 132-145 (2016).
15. Bulatov O. V., Elizarova T. G. Regularized shallow water equations for numerical simulation of flows with a moving shoreline. Comput. Math. and Math. Phys. 56, 661-679 (2016).
16. Потапов И. И., Снигур К. С. О решении уравнения Экснера для дна, имеющего сложную морфологию. Компьютерные исследования и моделирование 11 (3), 449-461 (2019).
17. Khrapov S. S., Khoperskov A. V. Application of graphics processing units for self-consistent modelling of shallow water dynamics and sediment transport. Lobachevskii Journal of Mathematics 41, 1475-1484 (2020).
18. Сухинов А. И., Чистяков А. Е., Сидорякина В. В. Проценко С. В., Атаян А. М. Локально-двумерные схемы расщепления для параллельного решения трехмерной задачи транспорта взвешенного вещества. Математическая физика и компьютерное моделирование 24 (2), 38-53 (2021).
19. Храпов С. С. Численное моделирование самосогласованной динамики поверхностных и грунтовых вод. Математическая физика и компьютерное моделирование 24 (3), 45-62 (2021).
20. Khrapov S., Pisarev A., Kobelev I., Zhumaliev A., Agafonnikova E., Losev A., Khoperskov A. The numerical simulation of shallow water: estimation of the roughness coefficient on the flood stage. Advances in Mechanical Engineering 2013 (id. 787016), 1-11 (2013).
21. Ландау Л. Д., Лифшиц Е. М. Теоретическая физика: гидродинамика. Т. 6. Москва, Наука (1986).
22. Exner F. M. Uber die Wechselwirkung zwischen Wasser und Geschiebe in Flussen Sitzungsber. Akad. Wiss. Wien, Math. Naturwiss. Kl. Abt. 2A 134, 165-180 (1925).
23. Liu X., Garcia M. Numerical investigation of seabed respone under waves with free-surface water flow. International Journal of Offshore and Polar Engineering 17 (2), 97-104 (2007).
24. Jayaratne R., Takayama Y., Shibayama T. Applicability of suspended sediment concentration formulae to large-scale beach morphological changes. Engineering Proceedings 1 (33), 57 (2012).
25. Зиновьев А. Т., Марусин К. В., Шибких А. А., Шлычков В. А., Затинацкий М. В. Математическое моделирование динамики течения и русловых деформаций на участке р. Обь у г. Барнаула. Ползуновский вестник 2, 204-209 (2006).
26. Harten A. High Resolution Schemes for Hyperbolic Conservation Laws. J. Comput. Phys. 49, 357-393 (1983).
27. Harten A., Lax P., van Leer B. On upstream differencing and Godunov type methods for hyperbolic conservation laws. SIAM review 25, 35-61 (1983).
28. Lu C. N., Qiu J.X., Wang R. Y. A numerical study for the performance of the WENO schemes based on different numerical fluxes for the shallow water equations. Journal of Computational Mathematics 28 (6), 807-825 (2010).
29. Хиеу Ле М., Тхань Д. Н. Х., Прасат В. Б. С. Монотонные разностные схемы второго порядка аппроксимации на неравномерных сетках для двумерного квазилинейного параболического уравнения типа конвекции-диффузии. Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия 7 (65), вып. 2, 343-355 (2020). https://doi.org/10.21638/11701/spbu01.2020.216
30. Семенова Н. Н., Терлеев В. В., Сухорученко Г. И., Орлов Е. Е., Орлова Н.Е. Об одном методе численного решения системы параболических уравнений. Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия 3 (61), вып. 2, 230-240 (2016). https://doi.org/10.21638/11701/spbu01.2016.207
Статья поступила в редакцию 26 ноября 2021 г.;
доработана 7 июня 2022 г.; рекомендована к печати 17 ноября 2022 г.
Контактная информация:
Храпов Сергей Сергеевич — канд. физ.-мат. наук, доц.; khrapov@volsu.ru
Numerical modeling of hydrodynamic accidents: Erosion of dams and flooding of territories*
S. S. Khrapov
Volgograd State University, 100, pr. Universitetskii, Volgograd, 400062, Russian Federation
For citation: Khrapov S. S. Numerical modeling of hydrodynamic accidents: Erosion of dams and flooding of territories. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2023, vol. 10(68), issue 2, pp. 357-373. https://doi.org/10.21638/spbu01.2023.215 (In Russian)
A mathematical and numerical model of the joint dynamics of surface water and traction sediment is built, which takes into account the nonlinear dynamics of fluid and bottom deformation. The dynamics of surface waters is described by the equations of Saint-Venant, taking into account the spatially inhomogeneous distribution of the terrain. The transport of sediment loads is described by the nonlinear Exner equation, generalized to the case of a spatially inhomogeneous distribution of the parameters of the underlying surface. For the numerical integration of the Saint-Venant and Exner equations, a stable and well-tested CSPH-TVD method of the second order of accuracy is used, the parallel CUDA algorithm
*The authors thank the Ministry of Science and Higher Education of the Russian Federation for financial support (government task no. 0633-2020-0003).
of which is implemented as a software package "EcoGIS-Simulation" for high-performance computing on supercomputers with graphic coprocessors (GPU). Hydrodynamic modeling of the processes of erosion of the enclosing dam of a real hydraulic facility and flooding of adjacent territories has been carried out. The parameters of the penetration of the enclosing dam and flooding zones, formed as a result of the development of a hydraulic accident at the tailing dump, have been determined. Based on the results obtained, it was concluded that the proposed method for numerical modeling of the joint dynamics of surface water and traction sediment can be more universal.
Keywords: mathematical modeling, Saint-Venant and Exner equations, computational fluid dynamics, CSPH-TVD method, flood zones, shallow water, entrained sediment, dam breakthrough and erosion, tailing dump, parallel computing, CUDA algorithm.
References
1. Levi I.I. Dynamics of channel processes. Leningrad, Gosenergoizdat Publ. (1957). (In Russian)
2. Karaushev A. V. River hydraulics. Leningrad, Gosenergoizdat Publ. (1969). (In Russian)
3. Richards K. J., Taylor P. A. A numerical model of flow over sand waves in water of finite depth. Geophys. J. Int. 65, 103-128 (1981).
4. Van Rijn L. C. Sediment transport, Part I: Bed load transport. Journal of Hydraulic Engineering 110 (10), 1431-1456 (1984).
5. Petrov P. G. Motion of a bed load. J. Appl. Mech. Tech. Phys. 32 (5), 717-721 (1991).
6. Venditti J. G., Church M. A., Bennett S. J. Bed form initiation from a flat sand bed. J. Geophys. Res.: Earth Surface 110 (F01009), 1-19 (2005).
7. Alan G.D., Peter D.Th. Advances in the study of moving sediments and evolving seabeds. Surveys in Geophysics January 29 (1), 1-36 (2008).
8. Petrov A. G., Potapov I. I. Sediment transport under normal and tangential bottom stresses with the bottom slope taken into account. J. Appl. Mech. Tech. Phys. 55 (5), 812-817 (2014).
9. Temporary guidelines for the calculation of zones in case of a sudden breakthrough of the tailing dams. Belgorod, VIOGEM Publ. (1981). (In Russian)
10. Safety of hydraulic structures at industrial and energy facilities: Collection of documents. Moscow, ZAO STC PB Publ. (2010). (In Russian)
11. Khrapov S.S., Khoperskov A. V., Kuzmin N.M., Pisarev A. V., Kobelev I. A. Numerical scheme for modeling the dynamics of surface waters based on the combined SPH-TVD approach. Computational Methods and Programming: New Computational Technologies 12 (1), 282-297 (2011). (In Russian)
12. Bulatov O. V., Elizarova T.G. Regularized shallow water equations and an efficient method for numerical simulation of shallow water flows. Comput. Math. and Math. Phys. 51, 160-173 (2011).
13. Sukhinov A. I., Chistyakov A.E., Semenyakina A. A., Nikitina A. V. Parallel implementation of the tasks of transporting substances and restoring the bottom surface on the basis of schemes of increased order of accuracy. Computational methods and programming: new computational technologies 16 (2), 256-267 (2015). (In Russian)
14. Dyakonova T., Khoperskov A., Khrapov S. Numerical model of shallow water: the use of NVIDIA CUDA Graphics processors. Communications in Computer and Information Science 687, 132-145 (2016).
15. Bulatov O. V., Elizarova T. G. Regularized shallow water equations for numerical simulation of flows with a moving shoreline. Comput. Math. and Math. Phys. 56, 661-679 (2016).
16. Potapov 1.1., Snigur K. S. On the solution of the Exner equation for a bottom with a complex morphology. Computer research and modeling 11 (3), 449-461 (2019). (In Russian)
17. Khrapov S. S., Khoperskov A. V. Application of Graphics Processing Units for Self-Consistent Modelling of Shallow Water Dynamics and Sediment Transport. Lobachevskii Journal of Mathematics 41, 1475-1484 (2020).
18. Sukhinov A. I., Chistyakov A. E., Sidoryakina V. V. Protsenko S. V., Atayan A. M. Local-two-dimensional splitting schemes for parallel solution of the three-dimensional problem of suspended matter transport. Math. physics and computer simulation 24 (2), 38-53 (2021). (In Russian)
19. Khrapov S. S. Numerical modeling of self-consistent dynamics of surface and ground waters. Math. physics and computer simulation 24 (3), 45-62, (2021). (In Russian)
20. Khrapov S., Pisarev A., Kobelev I., Zhumaliev A., Agafonnikova E., Losev A., Khoperskov A. The numerical simulation of shallow water: estimation of the roughness coefficient on the flood stage. Advances in Mechanical Engineering 2013 (id. 787016), 1-11 (2013).
21. Landau L. D., Lifshitz E. M. Theoretical Physics: Hydrodynamics. Vol. 6. Moscow, Nauka Publ. (1986). (In Russian)
22. Exner F. M. Uber die Wechselwirkung zwischen Wasser und Geschiebe in Flussen Sitzungsber. Akad. Wiss. Wien, Math. Naturwiss. Kl. Abt. 2A 134, 165-180 (1925).
23. Liu X., Garcia M. Numerical investigation of seabed respone under waves with free-surface water flow. International Journal of Offshore and Polar Engineering 17 (2), 97-104 (2007).
24. Jayaratne R., Takayama Y., Shibayama T. Applicability of suspended sediment concentration formulae to large-scale beach morphological changes. Engineering Proceedings 1 (33), 57 (2012).
25. Zinoviev A. T., Marusin K. V., Shibkikh A. A., Shlychkov V. A., Zatinatskiy M. V. Mathematical modeling of flow dynamics and channel deformations in the Ob near the city of Barnaul. Polzunovsky Bulletin 2, 204-209 (2006). (In Russian)
26. Harten A. High Resolution Schemes for Hyperbolic Conservation Laws. J. Comput. Phys. 49, 357-393 (1983).
27. Harten A., Lax P., van Leer B. On upstream differencing and Godunov type methods for hyperbolic conservation laws. SIAM Review 25, 35-61 (1983).
28. Lu C. N., Qiu J.X., Wang R. Y. A numerical study for the performance of the WENO schemes based on different numerical fluxes for the shallow water equations. Journal of Computational Mathematics 28 (6), 807-825 (2010).
29. Hieu L. M., Thanh D. N., Surya Prasath V. B. Second order monotone difference schemes with approximation on non-uniform grids for two-dimensional quasilinear parabolic convection-diffusion equations. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy 7 (65), iss. 2, 343-355 (2020). https://doi.org/10.21638/11701/spbu01.2020.216 (In Russian) [Eng. transl.: Vestnik St Petersburg University. Mathematics 53, iss. 2, 232-240 (2020). https://doi.org/10.1134/S1063454120020107].
30. Semenova N.N., Terleev V.V., Suhoruchenko G.I., Orlova E.E., Orlova N. E. On one method for the numerical solution of a system of parabolic equations. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy 3 (61), iss. 2, 230-240 (2016). https://doi.org/10.21638/11701/spbu01.2016.207 (In Russian) [Eng. transl.: Vestnik St Petersburg University. Mathematics 49, iss. 2, 138-146 (2016). https://doi.org/10.3103/S1063454116020138].
Received: November 26, 2021 Revised: June 7, 2022 Accepted: November 17, 2022
Author's information:
Sergey S. Khrapov — khrapov@volsu.ru
