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

Дьяконова Татьяна Андреевна аспирант, кафедра информационных систем и компьютерного моделирования, Волгоградский государственный университет 400062, Россия, г. Волгоград, пр. Университетский, 100 Dyakonova Tatyana post-graduate student of the Department of Information System and Computer Simulation at Volgograd State University 400062, Russia, Volgograd, pr. Universitetskii, 100 dyakonova_10@mail.ru Хоперсков Александр Валентинович доктор физико-математических наук профессор, кафедра информационных систем и компьютерного моделирования, Волгоградский государственный университет 400062, Россия, г. Волгоград, пр. Университетский, 100 Khoperskov Alexander Doctor of Physics and Mathematics professor of the Department of Information Systems and Computer Simulation at Volgograd State University 400062, Russia, Volgograd, pr. Universitetskii, 100 khoperskov@rambler.ru

Храпов Сергей Сергеевич кандидат физико-математических наук доцент, кафедра информационных систем и компьютерного моделирования, Волгоградский государственный университет 400062, Россия, г. Волгоград, пр. Университетский, 100 Khrapov Sergei PhD in Physics and Mathematics associate professor of the Department of Information Systems and Computer Simulation at Volgograd State University 400062, Russia, Volgograd, pr. Universitetskii, 100 xss-ip@mail.ru

Работа выполнена при финансовой поддержке РФФИ (гранты 16-07-01037, 15-45-02655)

Методы компьютерного моделирования динамики поверхностных вод в приближении мелкой воды позволяют решать широкий круг задач, связанный с прогнозом и управлением гидрологического режима на заданной территории с учетом разливов рек и озер ^[1], осадков ^[2-4], взаимодействия руслового и пойменного потоков ^[5], аварийных ситуаций на гидросооружениях ^{[6, 7]}, возникновением и динамикой цунами ^[8-10], формированием погоды вблизи крупных водоемов, морскими ^[11] и океаническими движениями ^[12], проведением экспертизы инженерных сооружений ^[13]. На основе многокомпонентных гидродинамических моделей мелкой воды с учетом динамики примесей (песка, ила, аллювиальными отложениями, транспорта донных наносов) возможно исследование изменений рельефа дна ^[14-16], меандрирования русел ^[17].

Плотины на реках, с одной стороны предназначены для регулирования водного стока, с другой представляют собой потенциальную угрозу территориям, расположенным вниз по течению, поэтому задача прогноза последствий затопления в случае чрезвычайных (аварийных) ситуаций на гидросооружениях является актуальной, и жизненно необходимой при разработке планов застройки и эвакуации населения при возникновении ЧС ^{[18, 19]}.

В данной работе описана математическая модель с учетом различных физических факторов и ее численная реализация для проведения вычислительных экспериментов на ЭВМ с графическими процессорами. Представлены результаты моделирования динамики затопления территории северной части Волго-Ахтубинской поймы вследствие аварии на плотине Волжской ГЭС.

Динамика поверхностных вод будем описывать в рамках модели мелкой воды (рис. 1), используя уравнения Сен-Венана ^[6], которые получаются после усреднения исходных трехмерных уравнений гидродинамики по вертикальной координате.

Рис. 1. Схема, поясняющая обозначения математической модели

Перечислим физические факторы, учитывающие в модели: 1) неоднородность рельефа дна `b(x,y)`); 2) силу тяжести `g`; 3) сила придонного трения; 4) взаимодействие между поверхностным потоком воды и приземным слоем атмосферы; 5) наличие осадков; 6) инфильтрация воды в подпочвенный слой; 7) испарение воды ^[20]. Уравнения Сен-Венана записываются относительно толщины слоя воды `H(x,y,t)` , горизонтальных компонент скорости, усредненных по вертикальной координате `u_x(x,y,t), u_y(x,y,t)`:

`(delH)/(delt)+(delHu_(x))/(delx)+(delH_(y))/(dely)=q,`

` `

(1)` `

`(delu_(x))/(delt)+u_(x)(delu_(x))/(delx)+u_(y)(delu_(x))/(dely)=-g(deleta)/(delx)+f_(x),`

(2)

`(delu_(y))/(delt)+u_(x)(delu_(y))/(delx)+u_(y)(delu_(y))/(dely)=-g(deleta)/(dely)+f_(y),`

(3)

где `eta=H+b` – определяет свободную поверхность воды, `f_(x), f_(y)` – соответствующие горизонтальные компоненты вектора `barf` сил, действующих на вертикально усредненный слой жидкости ^[20]. Для вектора `barf`можно записать:

`bar(f) = barf^"frict"+barf^"Cor"+barf^"wind"+barf^"visc",`

(4)

где `barf^"frict"=-lambda/2baru|baru|` – сила придонного трения,`lambda=(2gn_M^2)/(H^(4/3))` – коэффициент гидравлического сопротивления,`n_M` – коэффициент шероховатости по Маннингу ^{[16, 21-23]}, `barf^"Cor"=2[baruxxbbOmega]` - сила Кориолиса, `bbOmega` – угловая скорость вращения Земли,`barf^"wind"=C_arho_a/(rhoH)(barw-baru)|barw-baru|` – сила, обусловленная скоростью ветра `barw` (обычно на высоте 10 м от поверхности), `C_a` – параметр, характеризующий состояние водной поверхности (величина `C_a` характеризует своего рода «шероховатость» водной поверхности из-за наличия волн ^[24]), `rho_a` и `rho` плотность воздуха и воды соответственно,`bar(f)^"visc" = nu((del^2u_x)/(delx^2)+(del^2u_y)/(dely^2))` – сила внутреннего (вязкого) трения,`nu` – коэффициент кинематической турбулентной вязкости (`nu<=100` м/c²). Величина `n_M` является феноменологическим параметром, различаясь для конкретных водоемов в достаточно широких пределах `n_M~~0.01-:0.1` ^[25]. Специальные исследования для русла Волги ниже по течению Волжской плотины привели к значениям `n_M~~0.02-:0.07` в зависимости от уровня воды ^[26]. В общем случае величина `C_a` зависит от скорости ветра ^[27], однако во многих случаях можно пренебречь этой зависимостью, ограничившись `C_a=0.001` ^{[24, 28, 29]}. Поскольку во время затопления основные потоки воды распространяются вне пределов русел, то важным является учет пространственной неоднородности коэффициента Маннинга `n_M(x,y)` . Сходная ситуация возникает для задачи затопления в результате дождевых осадков ^[3].

В уравнении непрерывности для тонкого слоя (1) присутствует параметр `q` , который характеризует скорость притока/оттока воды за счет действующих источников/стоков. Для `q` в общем случае имеем:

`q=q^"(s)"+q^"(r)"-q^"(inf)"-q^"(ev)",`

(5)

где `q^"(s)"=(delQ)/(delS)` – приток воды через плотину, `Q(t)` – гидрограф попуска или зависимость скорости сброса воды через створ плотины (м³/с) от времени ^[25] (объем воды, попадающий в русло Волги за 1 секунду), `S` – площадь источника сброса, `q^"(r)"(t)` – приток воды, обусловленный осадками, `q^"(inf)"(x,y,t)` и `q^"(ev)"(t)` – отток воды за счет инфильтрации в почву и испарения соответственно. Указанные источники/стоки зависят от температуры воздуха, температуры воды, влажности воздуха, облачности, глубины промерзания почвы, насыщенности влагой грунта и его температуры.

Использование нелинейной модели инфильтрации воды в почву позволяет адекватно моделировать процесс впитывания воды ^[30]. Величину `q^"(inf)"` , входящую в соотношение (5), можно определить следующим образом:` `

`{(q^"(inf)"=q_"*"^"(inf)"(1-alpha)H/H_"*"),((delalpha)/(delt)=(q^"(inf)")/(psiH_l)-alpha/tau_l):}`

(6)

где `q_"*"^"(inf)"` – скорость впитывания воды толщиной слоя `H_"*"` в сухую почву, `alpha=V_w/V_l` – коэффициент влагонасыщенности почвы `(0<alpha<1)`, `V_w` – объем воды содержащийся в почвенном покрове с объемом `V_l` , `psi`– пористость почвы, `tau_l` – характерное время, за которое происходит осушение почвы за счет испарения и инфильтрации воды в подпочвенный слой с малой пористостью (глина в условиях Волго-Ахтубинской поймы).

Учитывая равнинный характер Волго-Ахтубинской поймы, перепад высот для которой лежит в пределах порядка 10 м, численное решение системы (1) – (3) можно строить на декартовой сетке с числом ячеек `N_(x)xxN_y` . При больших перепадах высот в горной местности выбор формы ячеек может заметно влиять на эффективность численной модели ^[3]. Размер ячейки `Deltax xx Deltay` ` ` ` `определяет погрешность расчета и выбирается исходя из качества рельефа `b(x,y)` , которую будем называть также функцией рельефа. Характерной особенностью рассматриваемой задачи является сложный характер рельефа, имеющий многочисленные нерегулярные и даже разрывные участки. Это предъявляет особые требования к численным алгоритмам интегрирования уравнений мелкой воды.

Основу алгоритма образует оригинальный численный метод, основанный на совместном использовании Эйлерова и Лагранжева подходов (TVD = Total-Variation-Diminishing и SPH = Smooth Particle Hydrodynamics соответственно) ^[31], важной особенностью которого является устойчивый сквозной счет для нестационарной границы типа «вода – сухое дно». Подробное описание данной численной схемы опубликовано в работах ^{[26, 31-33]}. Смешанные Лагранжево-Эйлеровые методы активно развиваются для решения различных задач ^{[34, 35]}. Реализованный численный алгоритм обеспечивает второй порядок точности интегрирования по времени и второй порядок точности по пространственным координатам. Численная схема является консервативной и хорошо сбалансированной.

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

Рис. 2. Иерархическая система сеток разных масштабов для Волго-Ахтубинской поймы

Ключевым моментом в моделировании и прогнозировании ЧС является своевременное получение информации для дальнейших действий по предотвращению или устранению последствий. Это требует использовать параллельные технологии для увеличении скорости вычислений, так, чтобы время проведения вычислительного эксперимента было меньше характерного времени рассчитываемого события. Наиболее перспективной представляется технология CUDA, позволяющая проводить вычисления на графических процессорах ускорителей. Проблемы построения, повышения эффективности и использования производительных вычислительных ресурсов невысокой стоимости и суперкомпьютеров активно обсуждаются в литературе ^{[4, 36-40]}. В данной работе в качестве вычислительных машин использовались рабочие станции с графическими ускорителями Tesla K40 и Tesla K80.

Использование подхода иерархической системы сеток (ИСС) хорошо адаптируется под программирование для GPU, поскольку в соответствии с архитектурой CUDA выполнение вычислительных ядер организовано как сеть потоковых блоков (рис. 3), которые являются аналогом блоков ИСС.

Рис. 3. Иерархия нитей на GPU

В настоящее время все большую популярность набирают гетерогенные вычислительные системы ^[41-43], которые состоят из многоядерных процессоров и массивно-параллельных ускорителей. Использование таких высокопроизводительных систем позволяет моделировать сложные нестационарные многомерные физические процессы, требующие больших затрат ресурсов. Поэтому при создании приложений для таких систем используются схемы, основанные на сочетании различных параллельных технологий. Для суперкомпьютеров с различным числом GPU (обозначим количество графических ускорителей величиной M) применяется двухуровневая параллельная модель OpenMP-CUDA, сочетающая в себе принципиально разные типы параллелизма (рис. 4a). Поскольку при загрузке GPU в кластере используется отдельный поток, для таких машин удобно использовать технологию OpenMP совместно с CUDA. Для обмена данными между различными GPU без копирования на CPU следует использовать доступ к памяти Direct Access, влияющий на повышение производительности.

Рис. 4. Гибридные схемы распараллеливания а) OpenMP-CUDA; б) MPI-OpenMP-CUDA

Для распараллеливания программ на гетерогенных системах, содержащих несколько CPU и GPU, используется трехуровневая модель, представленная на рис. 4б, которая сочетает в себе технологии MPI, OpenMP и CUDA. Применение гибридной схемы существенно ускоряет работу программ.

Рис. 5. Диаграмма деятельности расчетного модуля

Численное моделирование прорыва разрушения плотины проводилось на CUDA-версии расчетного модуля, диаграмма деятельности которой представлена на рис. 5. На диаграмме были использованы следующие обозначения:

- K_1 – расчет гидродинамических сил в момент времени t_n и шага по времени dt;

- K_2 – расчет интегральных характеристик и положения частиц в момент времени t_(n+1/2);

- K_3 – расчет гидродинамических сил в момент времени t_(n+1/2);

- K_4 – расчет интегральных характеристик и положения частиц в момент времени t_(n+1);

- K_5 – вычисление потоков.

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

Рис. 6. Диаграмма деятельности версии OpenMP-CUDA

Диаграмма представлена для вычислительной станции, содержащей 4 графических ускорителя. Обозначения D_1, D_2, D_3, D_4 соответствуют номеру GPU.

Волжская ГЭС является одной из крупнейших станций большой мощности и располагается в нижнем течении р. Волги (рис. 7). К счастью, прорывы на крупных плотинах и дамбах происходят крайне редко, однако такая угроза реально существует. В качестве примеров можно указать на аварию на Саяно-Шушенской ГЭС (2009 г.), плотины Сент-Френсис высотой 59 м (1928 г.), на реке Рейран (1959 г.), Вайонт в итальянских Альпах (262 метра, 1963 г.), ГЭС Байньцяо в Китае (118 метров, 1975 г.), на реке Какве (1993 г.), на реке Цяньтан (2011 г.). Поэтому рассмотрение возможных последствий нестандартных ситуаций ниже по течению от плотины Волжской ГЭС имеет практический интерес ^[6], а разработанное программное обеспечение может являться основой для построения схемы эвакуации при необходимости в зависимости от внешних факторов.

Рис. 7. Северная часть Волго-Ахтубинской поймы

Моделирование проводилось для территории Волго-Ахтубинской поймы. Ключевым фактором, определяющим качество гидродинамического расчета, является точность цифровой модели рельефа `b(x,y)` , ее соответствие реальной местности ^[44]. Для построения сеточной функции `b_(ij)=b(x_i,y_j)` использовались три источника: данные дистанционного зондирования Земли (ДЗЗ); лоции дна; результаты собственных GPS-измерений.

В основе модели рельефа лежат спутниковые данные ASTER и SRTM, способные обеспечить разрешение до 20 м в плоскости Земли и до 1 м по вертикальной координате. Модель речного дна уточнялась по данным лоций. На основе GPS-измерений были добавлены протоки и ерики в пойме. На рис. 8 представлена актуализированная цифровая модель рельефа Волго-Ахтубинской поймы.

Рис. 8. Цифровая модель рельефа Волго-Ахтубинской поймы, используемая для моделирования (шкалы указаны в метрах)

В качестве начального распределения воды и скорости течения выбрано гидрологическое состояние ВАП соответствующее пику паводка 2011 г. ^[45]. Выбор такого состояния на максимуме весеннего попуска воды является наиболее неблагоприятным при исследовании последствий ЧС. Как видно из рис. 10, в первые минуты после прорыва формируется волна высотой 7–10 м в виде гидравлического скачка, распространяющася сначала со скоростью порядка 70 км/ч по руслу Волги, а затем из-за наличия уклона рельефа в сторону поймы, а также особенностей русла Волги происходит смещение направления основного потока в волне прорыва. Наиболее разрушительной волна прорыва оказывается для пойменного участка в окрестностях населенных пунктов Краснослободска и Средней Ахтубы. Глубина воды в момент прохода волны по данному участку составляет в среднем 7 м. Следует отметить, что данная территория наиболее плотно застроена и заселена. После выхода волны затопления на территорию поймы ее скорость заметно снижается и в среднем составляет 10–15 км/ч, а толщина слоя воды на затопленных участках в среднем равна 5 м.

Рис. 9. Динамика затопления ВАП при сильном аварийном сбросе (прорыве плотины) в различные моменты времени

Построенная компьютерная модель позволяет прогнозировать последствия развития аварийных или катастрофических событий в Волго-Ахтубинской пойме, а также проводить техническую экспертизу гидросооружений. Результаты моделирования могут использоваться различными службами для принятия оперативных решений при чрезвычайных ситуациях. Даже в случае предельно сильного уровня прохождения воды через створ Волжской плотины правобережье (урбанизированная территория, включающая Волгоград, Волжский, Южную промзону, Светлый Яр) не будет затронута затоплением за исключением некоторых низких участков. В отличие от всей поймы, которая вся способна оказать под водой. Отметим, что в работе ^[6] ранее рассматривалась задача затопления территории Волго-Ахтубинской поймы при полном разрушении плотины Волжской ГЭС и был сделан вывод о наличии опасности для южной части г. Волгограда. Однако проведенные здесь вычислительные эксперименты с использованием более реалистичных моделей цифрового рельефа местности и гидрологического сопротивления позволяют пересмотреть выводы работы ^[6].

Использование распараллеленного кода для графических процессоров на основе технологии CUDA позволяет сократить длительность вычислительного эксперимента по изучению динамики затопления поймы до нескольких часов при использовании процессоров типа Tesla K40.