CC BY
26
Б01: 10.15514/КРКА8-2021-33(4)-16
Калибровка к - £ модели турбулентности в
пакете OpenFOAM с помощью методов машинного обучения для моделирования потоков на склонах гор на основе эксперимента
Д.И. Романова, ОЯСЮ: 0000-0002-5771-4114 <romanovadi@gmail.com> Институт системного программирования им. В.П. Иванникова РАН,
109004, Россия, г. Москва, ул. А. Солженицына, д. 25 Московский государственный университет имени М.В. Ломоносова, 119991, Россия, Москва, Ленинские горы, д. 1
Аннотация. В работе проводится калибровка к - 8 модели турбулентности для потоков со свободной поверхностью в русле и на склоне. Для калибровки модели проводится эксперимент в наклонном лотке постоянного уклона прямоугольного сечения. В эксперименте с помощью трубки Пито измеряются значения давления в потоке на различном расстоянии от дна, с помощью которых так же получается профиль скорости потока. На основании данных эксперимента с использованием алгоритма оптимизации Нелдера-Мида производится калибровка к - 8 модели турбулентности. Откалиброванная модель турбулентности далее применяется для расчёта прорыва ледникового озера, вблизи ледника Малый Азау.
Ключевые слова: математическое моделирование; численное моделирование; снежная лавина; грязекаменный сель; склоновые потоки; ОрепРОАМ; МегРоат; многофазный поток; турбулентный поток; к-8 модель турбулентности; оптимизация; метод Нелдера-Мидаб эксперимент; Малый Азау
Для цитирования: Романова Д.И. Калибровка к - 8 модели турбулентности в пакете ОрепРОАМ с помощью методов машинного обучения для моделирования потоков на склонах гор на основе эксперимента. Труды ИСП РАН, том 33, вып. 4, 2021 г., стр. 227-240. DOI: 10Л5514ДSPRAS-2021-33(4)—16
Благодарности: Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 19-31-90105. Автор выражает благодарность сотрудникам НИИ Механики МГУ Иванову Олегу и Трифонову Владимиру, а также Коровиной Дарье, Гинзбург Нике, и Гинзбург Борису за предоставленные экспериментальные данные.
Experiment based calibration of k - £ turbulence model in OpenFOAM package for mountain slope flows using machine learning
techniques
D.I. Romanova, ORCID: 0000-0002-5771-4114 <romanovadi@gmail.com> Ivannikov Institute for System Programming of the Russian Academy of Sciences, 25, Alexander Solzhenitsyn st., Moscow, 109004, Russia. Lomonosov Moscow State University, GSP-1, Leninskie Gory, Moscow, 119991, Russia
Abstract. We calibrate the k - s turbulence model for free surface flows in the channel or on the slope. To calibrate the turbulence model, an experiment is carried out in an inclined rectangular research tray. In the experiment, the pressure values in the flow are measured at different distances from the bottom using a Pitot tube; after transforming data, the flow velocity profile is obtained. The k - s turbulence model is calibrated based on experimental data using the Nelder-Mead optimization algorithm. The calibrated turbulence model is then used to calculate the outburst of a lake near the glacier Maliy Azau on the Elbrus (Central Caucasus).
Keywords: mathematical modeling; numerical modeling; snow avalanche; mudflow; slope flow; OpenFOAM; interFoam; multiphase flow; turbulent flow; k-s turbulence model; Nelder-Mead optimization; free surface flow; Maliy Azau.
For citation: Romanova D.I. Experiment based calibration of k - s turbulence model in OpenFOAM package for mountain slope flows using machine learning techniques. Trudy ISP RAN/Proc. ISP RAS, vol. 33, issue 4, 2021, pp. 227-240 (in Russian). DOI: 10.15514/ISPRAS-2021-33(4)-16
Acknowledgements. The reported study was funded by RFBR, project number 19-31-90105. The author is grateful to the staff of the Research Institute of Mechanics of Moscow State University, Oleg Ivanov and Vladimir Trifonov, as well as Daria Korovina, Nika Ginzburg, and Boris Ginzburg for the experimental data provided.
1. Введение
Математическое моделирование динамики потоков на горных склонах позволяет прогнозировать опасные зоны, изучать размыв склонов или каналов потоками, а также рассчитывать напряжения, возникающие в защитных сооружениях и других объектах при ударе потока.
Большая часть используемого в настоящее время программного обеспечения для моделирования течений на горных склонах основывается на уравнениях мелкой воды (осреднённые по глубине уравнения механики сплошной среды). Такой подход дает информацию только о границах опасной зоны, толщине потока и средней скорости потока по глубине. Кроме того, осреднённые по глубине уравнения содержат коэффициенты, которые можно найти только при региональной калибровке. Трехмерное моделирование, основанное на полных уравнениях, а не на усредненных по глубине уравнениях, может предоставить полную информацию о распределении параметров в потоке. Это очень важно, особенно, в задаче взаимодействия потока с разными объектами.
Исторически первыми для исследования потоков на склонах гор использовались одномерные осреднённые по глубине модели в силу своей вычислительной простоты. Это такие работы как [2, 5, 11, 12, 14, 18, 34, 35, 39], датирующиеся начиная с 1967 года. Дальнейшим развитием осреднённых по глубине моделей было рассмотрение двумерного подхода [13, 25, 29, 40]. В настоящее время использование осреднённых по глубине уравнений можно видеть в работах [6-8, 17].
Несколько исследовательских групп в мире занимались моделированием потоков на склонах гор с использованием 3Б-моделирования [3, 4, 9, 15, 21, 22, 24, 30-34, 37, 41], среди них Ямагучи (Yuya Yamaguchi) и соавторы [41], они использовали метод конечных элементов на
основе подхода Петрова-Галёркина для моделирования склонового потока (снежная лавина) как двухфазного потока (снег и воздух), где снег был представлен бингамовской жидкостью. Кеничи Ода (Kenichi Oda) и соавторы [30] также использовали двухфазный подход для описания динамики снежной лавины, только в отличие от работы [41] система уравнений решалась конечно -разностным методом. В настоящей работе будет использоваться многофазный односкоростной подход для моделирования потоков на склонах гор с использованием метода конечного объёма [31]. Течение считается турбулентным, используется k - е модель турбулентности. Система уравнений решается методом конечного объёма.
Рис. 1. Схема экспериментального лотка [1] Fig. 1. Scheme of experimental chute [1]
0.175 0.150 0.125
f алоэ j=
□1
I 0.0 75 0.050 0.025
0.000
0 2 4 6 8 10
time (s)
Рис. 2. Графики скорости (синий) и глубины потока (красный), замеренные на расстоянии 11.1 метра от начала установки. Показаны экспериментальные значения из работы [1 ] (Experiment), вычисленные с использованием осреднённого по глубине подхода (Calculated d.-av.), вычисленные с использованием трёхмерного подхода при ламинарном режиме течения (Calculated lam), вычисленные
с использованием k - е модели турбулентности в трёхмерном подходе (Calculated turb. k-e) и вычисленные с использованием k - m SSTмодели турбулентности в трёхмерном подходе (Calculated
turb. k-w SST)
Fig. 2. Depth-averaged velocity (blue) and flow depth (red) graphs measured at a distance of 11.1 meters from the start of the installation. Shown are the experimental values from [1 ] (Experiment), calculated using depth-avaraged approach (Calculated d.-av.), calculated for laminar flow (Calculated lam), calculated using the k - е turbulence model (Calculated k-e) and calculated using the k - m SST turbulence model (Calculated
k-w)Scheme of experimental chute [1]
2. Актуальность
Потоки на склонах гор являются турбулентными потоками. Существующие модели турбулентности содержат в себе ряд констант, которые были откалиброваны при расчёте канонических течений таких, как обтекания потоком воздуха различных профилей, течения в трубах и каналах и др. [26, 27, 38]. Калибровки моделей турбулентности для расчёта
потоков неньютоновской жидкости со свободной поверхностью на склонах гор не проводилось.
На примере расчёта эксперимента по спуску водоснежного потока в экспериментальном лотке с двумя небольшими защитными дамбами и одной основной удерживающей дамбой (рис. 1), проведённого в университете Исландии [1,23], показана необходимость калибровки турбулентной модели с целью увеличения точности расчёта защитных сооружений (рис. 2 и табл. 1).
Данный эксперимент, поставленный в университете Исландии был спроектирован с целью воспроизвести взаимодействие водоснежного потока с различными препятствиями. Для проведения эксперимента использовалась чистая водопроводная вода, так как она является легко доступным материалом и схожа по плотности с насыщенным водоснежным потоком. В том числе вода позволяет воспроизвести отличия в динамике водоснежных потоков от сжимаемых потоков гранулированных сред с высоким числом Фруда. Эксперимент был спроектирован таким образом, чтобы число Фруда было близко к фактическому числу Фруда для водоснежных потоков: в диапазоне 3 < Fr < 5. Глубина потока, высота барьеров и другие геометрические характеристики эксперимента моделировались с коэффициентом 10. Эксперимент проводился в длинном деревянном желобе, который состоял из четырех секций с разными углами наклона (55° , 27° , 11.7° , 0° — слева направо), см. рис. 2. Самые нижние 5.05 м желоба были горизонтальными, а верхние его части были покрыты ковром для увеличения трения и снижения скорости потока и, следовательно, числа Фруда для получения целевого диапазона чисел Фруда. Защитные конструкции были выполнены из фанеры.
Табл. 1. Сравнение измеренных и рассчитанных параметров потока
Table. 1. Comparison o of measured and calculated flow parameters
Параметр сравнения Экспериментальные данные Данные численного расчёта
Ламинарный режим k - s модель турбулентно сти k - w SST модель турбулентно сти
Высота всплеска 1.3 м 2.28 м 1.64 м 2.4 м
на основной дамбе
Время взаимодействия потока с основной дамбой 1.25 с 1.5 с 1.3 с 1.3 с
Объём, 2.684 м3 2.629 м3 2.528 м3 2.655 м3
удержанный дамбой из 2.7 м3
Расчёт эксперимента был проведён с использованием многофазного односкоростного подхода, реализованного в пакете OpenFOAM в решателе interFoam. На рис. 2 видно, что расчёты с использованием турбулентных моделей лучше воспроизводят эксперимент нежели ламинарный режим течения, незначительное преимущество у k - s модели турбулентности. Однако по результатам сравнении параметров в таблице 1 вытекает, что оба турбулентных расчёта требуют доработки.
В данной работе с помощью методов машинного обучения производится калибровка k - s модели турбулентности для такого класса задач, как потоки со свободной поверхностью на склоне под действием силы тяжести.
3. Калибровочный эксперимент
Для калибровки коэффициентов к - £ модели турбулентности используется эксперимент, поставленный в НИИ Механики МГУ. В процессе эксперимента в лотке постоянного уклона (эксперимент проведён для нескольких разных углов наклона) производился спуск турбулентного потока жидкости. Использовался лоток следующей геометрии: длина — 1 м, ширина — 12 см, высота бортиков — 10 см.
ОоЛъ
_— "II
Рис. 4. Схема экспериментального лотка Fig. 4. Experiment chute scheme В расчётах использовалась область лотка между двумя точками замера скорости и глубины, находящимися на расстоянии 23 см и 82 см от верхнего края лотка. Схема эксперимента показана на рис. 4. Для проведения эксперимента используется водопроводная вода, как в эксперименте университета Исландии. Лоток выполнен полностью из акрилового стекла толщиной 4 мм. Вода подаётся из резервуара с возможностью регулировки высоты столба жидкости при подаче. При входе в исследуемый участок лотка установлен успокоитель потока. Эксперимент проводится в стационарном режиме с замкнутым жидкостным контуром. Стационарность обеспечивает насос перекачки жидкости из нижнего водосборника в верхний резервуар. Для замера профиля скорости потока используется трубка Пито, присоединённая к датчику давления.
Было проведено 3 серии экспериментов, в которых менялся угол наклона склона, начальная глубина потока и начальный профиль потока, как показано в табл. 2.
Табл. 2. Параметры расчётов Table. 2. Experiments' parameters
u0 среднее по глубине, м/с h0, мм в
1.63 4.20 25°
2.00 4.95 28°
1.78 3.45 33°
4. Математическая модель
Используется трёхмерный подход для моделирования потока. При этом подходе поток рассматривается как многофазное течение, осреднение по глубине не используется. Одна из фаз — воздух, другая — материал потока (снег или грязекаменная смесь), третьей фазой может быть материал подстилающей поверхности. В данной работе рассматривается двухфазная модель — поток-воздух. Используется к — е модель турбулентности, основанная на работах [27], [38].
В работе используется метод VOF (объём жидкости) для отслеживания границы свободной поверхности, который был предложен Хиртом (C.W. Hirt) и Николсом (B.D. Nichols) в 1981
году [20]. Данный метод не отслеживает границу явно, она задаётся как пороговое значение для объёмной доли фазы.
Выпишем систему уравнений для описания двухфазной модели течения, в которой каждая из фаз считается несжимаемой и обе фазы имеют единую скорость (1). Данная модель состоит из следующих уравнений: уравнения неразрывности для смеси, уравнения переноса объёмной доли одной из фаз, уравнения сохранения количества движения, уравнения для расчёта турбулентной кинетической энергии и уравнения диссипации турбулентной кинетической энергии, а также замыканий для тензора напряжений, являющегося функцией эффективной вязкости и тензора скоростей деформаций, эффективная вязкость является суммой молекулярной и турбулентной вязкостей, а турбулентная вязкость представлена функцией турбулентной кинетической энергии и её диссипации, плотность и вязкость смеси в модели вычисляются по принципу весового среднего.
V -и = 0,
да ~dt
+ V - (иа) = 0,
д(ри) __
д(рк) dt
д(ре) dt
д
+ V - (рии) = -Vp + V -Т + pf,
+ V - (рик) = V - (ßVk) - V - (—Vk + Рк) - pe,
\ак )
2
(1)
£ £2 /Ut \
+ V - (рие) = СЕ1Рк--pCE2 — +V - (—Ve),
т = 2neffs,s = 0.5[Fu + (Vu)T],ßeff = ß + ßt,ßt = pCß k2/e, p = p1a + p0(1 — a),ß = vp,v = v1a + v0(1 — a).
Здесь U — скорость смеси, горизонтальной чертой над буквами обозначается осреднение по Рейнольдсу; a — объёмная доля выбранной фазы; p — плотность смеси; р — давление; f — плотность массовых сил; k — плотность турбулентной кинетической энергии; е — диссипация турбулентной кинетической энергии; Рк — скорость производства турбулентной кинетической энергии средним течением, Се1, Се2, См, ак, аЕ — коэффициенты модели турбулентности.
Объёмная доля фазы принимает значения в диапазоне 0 < a < l. В случае, например, если a = 0 в ячейке, то она полностью заполнена фазой 0, или, если a = l в ячейке, то она полностью заполнена фазой 1.
Более детальное описание математической модели можно найти в книге Ферцигера (Joel Ferziger) и Перича (Milovan Peric) [16].
5. Вычислительная модель
Для решения гидродинамической задачи течения жидкости на склоне используется свободно распространяемый пакет с открытым исходным кодом OpenFOAM.
5.1. Расчётная область для эксперимента НИИ Механики МГУ
Для расчёта эксперимента НИИ Механики МГУ была рассчитана часть экспериментального лотка, находящаяся между двумя точками замера профиля скорости и глубины потока. Расчёт проводился для срединной полоски лотка толщиной 10 мм, где влияние боковых стенок незначительно. Первый измеренный профиль подавался на вход расчётной области, второй являлся объектом сравнения. Геометрия расчётной области представляет собой параллелепипед длинной 590 мм, шириной 10 мм, и глубиной 30 мм. Количество ячеек составило 590х10х60.
Были выделены следующие границы расчётной области: дно лотка, боковые стенки лотка, верхняя граница лотка, плоскость на входе в лоток, плоскость на выходе из лотка. Были заданы следующие граничные условия:
• дно лотка: для скорости задано условие прилипания, для объёмной доли фазы -нулевой градиент, для турбулентной кинетической энергии установлено условие нулевого градиента, для диссипации турбулентной кинетической энергии установлена константа
0.1 м2/с3;
• боковые стенки лотка: для всех рассматриваемых величин задано условие нулевого градиента;
• верхняя граница лотка: смешанное условие с заданием атмосферного давления, и условия отсутствия притока среды через данную границу, отток происходит по принципу нулевого градиента;
• входная плоскость: заданы фиксированные значения объёмной доли фазы и профиля скорости потока, для других величин установлен свободный вход с условием нулевого градиента, позволяющий с течением времени установиться исследуемым параметрам потока;
• выходная плоскость: для всех величин установлено условие нулевого градиента. Начальные условия в задаче таковы, что объём полностью заполнен неподвижным воздухом и подаётся жидкость через входную плоскость, спустя время поток устанавливается и снимаются замеры на выходной плоскости для сравнения с экспериментальными данными. Поток считается установившимся спустя 5 секунд.
5.2. Расчётная область для эксперимента, проведённого в университете Исландии
Для расчёта эксперимента, поставленного в университете Исландии использовалась двумерная расчётная область. Были выделены следующие границы:
• дно лотка, включающее в себя защитные сооружения, границу окончания экспериментальной установки и борта резервуара с жидкостью, для которых заданы аналогичные с экспериментом НИИ Механики МГУ граничные условия;
• верхняя граница экспериментальной установки, с соответствующими граничными условиями, описанными в предыдущем разделе;
• боковые границы расчётной области, заданные как нерасчётные для реализации двумерной постановки.
Всего расчётная область состоит из 257 тысяч ячеек. Начальные условия заданы таким образом, что вода, глубиной 0.9м располагается в стартовом резервуаре и удерживается виртуальной границей, исчезающий с началом расчёта.
6. Алгоритм оптимизации
Необходимость и эффективность калибровки турбулентных моделей продемонстрированы в таких работах, как [19,28]. Работа учёных Сержа Гийяса (Serge Guillas), Нины Гловер (Nina Glover) и Лиоры Малки-Эпштейн (Liora Malki-Epshtein) [19] демонстрирует эффективность применения калибровки с помощью байесовских методов для коэффициентов k - е модели турбулентности при моделировании воздушных потоков в городе. В работе американских учёных Джулии Линг (Julia Ling), Эндрю Курзавски (Andrew Kurzawski) и Джереми Темплтона (Jeremy Templeton) [28] авторы используют глубокие нейронные сети для получения замыкающего члена для расчёта напряжений Рейнольдса, учитывающего анизотропию турбулентности потока жидкости в канале. В применении к задаче моделирования динамики снежной лавины или грязе-каменного селя точность моделирования динамики турбулентного потока имеет сильное влияние на расчёт защитного
сооружения. Как было показано в табл. 1, различие экспериментальных и расчётных параметров может достигать 25%.
Рис. 3. Архитектура предлагаемого алгоритма для калибровки коэффициентов k - е модели
турбулентности
Fig. 3. The architecture of the proposed algorithm for calibrating the k - е turbulence model coefficients k - e модель турбулентности содержит ряд констант: C = 0.09; Ci = 1.44; C2 = 1.92; C3 = 0.0; Gk = 1.0; Ge = 1.3. Использование данной турбулентной модели для расчёта потоков неньютоновских жидкостей на склонах требует калибровки коэффициентов модели. Для организации процесса калибровки используются алгоритмы машинного обучения. Калибровка турбулентной модели проводилась с использованием следующего алгоритма, основанного на обучении с подкреплением, показанного на рис. 3:
1. Обучение на базе ряда расчётов, проведённых с использованием RANS модели с k-e моделью турбулентности с различными значениями констант;
2. Получение новых значений коэффициентов турбулентной модели посредством алгоритма оптимизации;
3. Расчёт гидродинамики потока с использованием турбулентной модели с коэффициентами, полученными с помощью алгоритма оптимизации;
4. Дообучение алгоритма с использованием полученных данных расчёта гидродинамики потока.
При оптимизации коэффициентов турбулентной модели минимизируется корень из среднеквадратичного отклонения вычисленного профиля скорости потока на выходной плоскости от экспериментального профиля:
M
£z=0 {fexp О) - vk (z)}
(2)
h
где И — глубина потока на выходной плоскости.
Для минимизации используется метод Нелдера-Мида (метод безусловной оптимизации функции от нескольких переменных, не использующий градиентов), реализованный в библиотеке SciPy.
Для реализации приведённого выше алгоритма было разработан модуль, который позволил в автоматическом режиме связать OpenFOAM расчёты с алгоритмом оптимизации. В результате в процессе оптимизации генерировались обновлённые значения коэффициентов турбулентной модели, которые передавались в утилиту препроцессинга, готовящую расчётные кейсы. Далее запускались расчёты OpenFOAM кейсов в параллельном режиме, по результатам расчётов рассчитывалось значение loss-функции и передавалось в оптимизатор.
7. Результаты калибровки
Для оптимизации коэффициентов к - е модели турбулентности было проведено 180 расчётов калибровочного эксперимента (рис. 5), что заняло 213 часов работы на 24 ядрах. Начальные значения коэффициентов были заданы следующими:
Си = 0.09; С1 = 1.44; С2 = 1.92; С3 = 0.0; ок = 1.0; аЕ = 1.3.
(3)
После калибровки значения коэффициентов стали следующими:
См = 0.0902; С1 = 1.5419; С2 = 1.92; С3 = 0.0001; ак = 1.0215; ае = 1.3334. (4) Были получены профили скорости на выходной плоскости для различных углов наклона лотка, как показано на рис. 5
Калибровка привела к минимизации функции потерь (5), показанной в табл. 3
— 25° Experiment — 25" Calculated k-e (initial) —- 25* Calculated k-e {tuned] [ — 28° Expemnent — 28* Calculated k-e (initial| — 33* Expriment — 33* Calculated k-e (initial]
| |
j ' ' i
I / П ' /
j if
чГ /t
Рис. 5. Сравнение экспериментального профиля скорости с расчётным при использовании стандартных значений коэффициентов k - е модели турбулентности и откалиброванных значений
коэффициентов в лотках различного угла наклона к горизонту Fig. 5. Comparison of the experimental velocity profile with the calculated one using the standard values of the k - е turbulence model coefficients and the calculated velocity profile with calibrated values of the coefficients for different slope inclination angles (to the horizon)
Табл. 3. Минимизация функции потерь
Угол наклона лотка Начальное значение функции потерь Минимизированное значение функции потерь
25° 0.242 0.218
28° 0.377 0.389
33° 0.236 0.168
8. Верификация на эксперименте университета Исландии
Был произведён расчёт эксперимента университета Исландии, показанного на рис. 1, с использованием откалиброванных коэффициентов к - е модели турбулентности. Были получены результаты, показанные на рис. 6 и в табл. 4.
На рис. 6 показано, что для глубины потока уменьшилось различие экспериментальных данных от расчётных после калибровки коэффициентов турбулентной модели. Для параметров, перечисленных в таблице 4 можем видеть уменьшение расхождения для высоты первого всплеска на основной дамбе на 41%, для объёма, удержанного дамбой — на 60%. Из вышесказанного можно сделать вывод, что калибровка коэффициентов турбулентной модели для исследуемого класса задач увеличила точность расчётов. Однако, не будет лишней дальнейшая калибровка турбулентной модели для учёта анизотропии турбулентности, которая имеет существенное значение при исследовании потоков со свободной поверхностью на склонах.
0.175 0.153 0.125
f 0.100
J= □1
0.0 75 0.050 0.025
0.000
О 2 4 6 8 10
time (s)
Рис. 6. Графики скорости (синий) и глубины потока (красный), замеренные на расстоянии 11.1 метра от начала установки. Показаны экспериментальные значения из [1] (Experiment), вычисленные с использованием k - е модели турбулентности (Calculated k-e (initial)) с исходными значениями коэффициентов и вычисленные с использованием k - е модели турбулентности (Calculated k-e (tuned)) с откалиброванными коэффициентами Fig. 6. Depth-averaged velocity (blue) and flow depth (red) graphs measured at a distance of 11.1 meters from the start of the installation. Shown are the experimental values from [1] (Experiment), calculated using the k - е turbulence model (Calculated k-e) with initial coefficients and calculated using the k - е turbulence model (Calculated k-e tuned) with calibrated coefficients
Табл. 4. Сравнение измеренных и рассчитанных параметров потока Table. 4. Loss function minimization
Параметр сравнения Экспериментальные данные k - s модель турбулентности с исходными коэффициентами k - s модель турбулентности с откалиброванными коэффициентами
Высота первого всплеска на основной дамбе 1.3 м 1.64 м 1.5 м
Время с начала взаимодействия потока с основной дамбой, до окончания переливания потока через дамбу 1.25 с 1.3 с 1.3 с
Объём, удержанный дамбой из 2.7 м3 2.684 м3 2.528 м3 2.622 м3
9. Расчёт возможного прорыва озёр вблизи ледника Малый Азау
С помощью полученной откалиброванной модели турбулентности произведём расчёт возможного прорыва ледникового горного озера вблизи ледника Малый Азау [10]. Группа озёр, расположенная около ледника Малый Азау на южном склоне Эльбруса, широко известна и ежегодно посещается туристами и альпинистами. В период 1957- 2021 гг. зафиксировано два случая прорывов озёр: озеро западное в 1978 г. и озеро восточное в 2011 г. Первый прорыв связан с процессами оползневых деформаций моренного массива,
слагающего часть озёрной котловины, второй — с деградацией ледяной плотины озера и переливом воды поверх неё. Существующие в настоящее время моренные плотины озёр и участки поверхностного стока воды из озёр находятся в стабильном состоянии, поэтому угроза прорыва озёр отсутствует. Потенциальная угроза прорыва остаётся в связи с высокой сейсмичностью и возможной вулканической деятельностью в этом районе. Разработанное ПО позволяет получить точные данные о пути движения паводка, его профиле скорости, распределении давления в потоке, напряжении на дне потока в различных точках.
Рис. 7. Карта рельефа исследуемой местности (слева) [10]; трёхмерная расчётная область в начальный момент времени с начальным расположением озера (центр); картина потока в момент
времени 100 минут от начала прорыва (справа) Fig. 7. Maliy Azau glacial lake outburst flood map (left) [10]; three-dimensional computational domain at the initial time moment with the initial location of the lake (center); outburst flood at 100 minutes from the
start (right)
Рис. 8. Графики расхода жидкости (красный) и объёма жидкости слившегося с озера (синий) Fig. 8. Graphs of the flow rate (red) and the volume of liquid poured out of the lake (blue) В результате математического моделирования прорыва озера глубины 2 метра были получены картины течения, как показано на рис. 7. Слева на рис. 7 показана карта рельефа
исследуемой области, зелёным цветом схематически отмечен контур расчётной области; в центре рис. 7 показана расчётная область в начальный момент времени, фиолетовым цветом отмечена вода в озере; справа на рис. 7 показано распределение скорости потока спустя 100 минут после прорыва.
На рис. 8 показано, что прорыв озера продолжается более 200 минут, максимальный расход жидкости составляет порядка 760 м3 /мин, средний — 58 м3 /мин. Также на рис. 8 мы видим всплеск расхода жидкости сразу после прорыва, характерный для такого типа течений и дальнейшее выравнивание потока.
10. Результаты и выводы
В процессе работы был проведён эксперимент, по результатам которого была откалибрована k —£ модель турбулентности, произведён расчёт возможного прорыва ледникового горного озера вблизи ледника Малый Азау (Кавказ).
• Было получено уменьшение расхождения на обучающем эксперименте до 29%.
• При оптимизации коэффициентов турбулентной модели использовались три эксперимента, что позволило избежать переобучения алгоритма.
• При использовании турбулентной модели с откалиброванными коэффициентами при моделировании тестового эксперимента Университета Исландии получено уменьшение расхождения расчётных характеристик с экспериментальными.
• При использовании турбулентной модели с откалиброванными коэффициентами при моделировании тестового эксперимента Университета Исландии получено уменьшение расхождения расчётного объёма, удержанного комбинацией защитных сооружений и экспериментального на 60%.
• С использованием откалиброванной модели турбулентности было проведено моделирование возможного прорыва ледникового горного озера вблизи ледника Малый Азау (Кавказ).
Список литературы / References
[1]. Katrin Helga Agustsdottir. The design of slushflow barriers: Laboratory experiments. PhD thesis, Faculty of Industrial Engineering, Mechanical Engineering, and Computer Science, University of Iceland, 2019.
[2]. Бахвалов Н.С., Эглит М.Э. Исследование решений уравнения движения снежных лавин. Материалы гляциологических исследований. Хроника обсуждений, вып. 16, 1970 г., стр. 31-38 / N.S. Bakhvalov, M.E. Eglit. Investigation of the solutions to snow avalanche movement equations. Glaciological data. GD-16, 1984, pp. 117-128.
[3]. Eloise Bovet, Bernardino Chiaia, Luigi Preziosi. A new model for snow avalanche dynamics based on bingham fluids. Meccanica, vol. 45, 2010, pp. 753-765.
[4]. W. Brandstatter, F. Hagen et al. Dreidimensionalle simulation von staublawinen unter berucksichtigung realer gelandeformen. Zeitschrift der Wildbach- und Lawinebnverbauung Österreichs, 120, 1992, pp. 107137 (in German).
[5]. A. Briukhanov, S. Grigorian et al. On some new approaches to the dynamics of snow avalanches. In Proc. of the Conference on Physics of Snow and Ice, 1966, pp, 1223-1241.
[6]. Yves Bühler, Marc Christen et al. Sensitivity of snow avalanche simulations to digital elevation model quality and resolution. Annals of Glaciology, vol. 52, 2011, pp. 72-80.
[7]. M. Christen, J. Kowalski et al. Numerical simulation of dense snow avalanches in three-dimensional terrain. Cold Regions Science and Technology, vol. 63, issue 1, 2010, pp. 1-14.
[8]. Marc Christen, Perry Bartelt et al. Calculation of Dense Snow Avalanches in Three-Dimensional Terrain with the Numerical Simulation Program Ramm. In Proc. of the Whistler 2008 International Snow Science Workshop, 2008, pp. 709-716.
[9]. J. Dent, T. Lang. Modeling of snow flow. Glaciology, vol. 26, 1980, pp. 131-140.
[10]. M.D Dokukin, A Khatkutov. Lakes near the glacier maliy azau on the elbrus (central caucasus): dynamics and outbursts. Ice and Snow, vol. 56, 2016, pp. 472-479.
[11]. М.Э. Эглит. Теоретические подходы к расчету движения снежных лавин. В сборнике Итоги науки. Гидрология суши, гляциология. М., ВИНИТИ, 1968, стр. 60-98 / M. Eglit. Theoretical approaches to the calculation of the motion of snow avalanches. Glaciological data. GD-16, 1984, pp. 63-118.
[12]. М.Э. Эглит. Вычисление параметров лавин в зонах торможения и остановки. Материалы гляциологических исследований, вып. 43, 1982 г., стр. 35-39 / M. Eglit. Calculation of the parameters of avalanches in the runout zone. Materialy Glyatsiologicheskikh Issledovaniy, issue 53, 1982, pp. 35-39 (in Russian).
[13]. М.Э. Эглит. Неустановившиеся движения в руслах и на склонах. Изд-во МГУ, 1986 г., 96 стр. / M. Eglit. Unsteady motions in channels and on slopes. Moscow State University, 1986, 96 p. (in Russian).
[14]. Margarita Eglit. Some mathematical models of snow avalanches. In Advances in the Mechanics and the Flow of Granular Materials, vol. 2, 1983, pp. 557-588.
[15]. M.E. Egli, A.E. Yakubenko. Numerical modeling of slope flows entraining bottom material. Cold Regions Science and Technology, vol. 108, 2014, pp. 139-148.
[16]. Joel Ferzige , Milovan Peric. Computational Methods for Fluid Dynamics, 3rd ed. Springer, 2001, 426 p.
[17]. Jan-Thomas Fischer, Julia Kowalski et al. Dynamic avalanche modeling in natural terrain. In Proc. of the International Snow Science Workshop, 2009, pp. 448-453.
[18]. С.С. Григорян, М.Э. Эглит, Ю.Л. Якимов. Новая математическая постановка задачи о движении лавины и решение этой задачи. Труды Высокогорного геофизического института, no. 12, 1967 г., стр. 104-113 / S.S. Grigorian, M.E. Eglit, Y.L. Iakimov. A new formulation and solution of the problem of snow avalanche motion. Trudy Vycokogornogo Geofizicheskogo Instituta, no. 12, 1967, pp. 104-113 (in Russian).
[19]. Serge Guillas, Nina Glover, Liora Malki-Epshtein. Bayesian calibration of the constants of the k - s turbulence model for a cfd model of street canyon flow. Computer Methods in Applied Mechanics and Engineering, vol. 279, 2014, pp. 536-553.
[20]. C.W. Hirt, B.D. Nichols. Volume of fluid (vof) method for the dynamics of free boundaries. Journal of Computational Physics, vol. 39, issue 1, 1981, pp. 201-225.
[21]. Dieter Issler. Modelling of snow entrainment and deposition in powder-snow avalanches. Annals of Glaciology, vol. 26, 1998, pp. 253-258.
[22]. Dieter Issler, Manuel Pastor Pérez. Interplay of entrainment and rheology in snow avalanches: a numerical study. Annals of Glaciology, vol. 52, issue 58, 2011, pp. 143-147.
[23]. Rebecca Anne Jones. The Design of Slushflow Barriers: CFD Simulations. PhD thesis, Faculty of Industrial Engineering, Mechanical Engineering, and Computer Science, University of Iceland, 2019.
[24]. M.A. Kern, F. Tiefenbacher, J.N. McElwaine. The rheology of snow in large chute flows. Cold Regions Science and Technology, vol. 39, issue 2-3, 2004, pp. 181-192.
[25]. A.G. Kulikovskii, M.E. Eglit. Two-dimensional problem of the motion of a snow avalanche along a slope with smoothly changing properties. Journal of Applied Mathematics and Mechanics, vol. 37, issue 5, 1973, pp. 792-803.
[26]. A. Morse, B.E. Launder et al. Prediction of free shearflows - a comparison of the performance of six turbulence models. In Proc. of the of NASA Conference on Free Shear Flows, 1972, pp. 361-426.
[27]. Brian Launder, D.B. Spalding. The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering, vol. 3, issue 2, 1974, pp. 269-289.
[28]. Julia Ling, Andrew Kurzawski, Jeremy Templeton. Reynolds averaged turbulence modelling using deep neural networks with embedded invariance. Journal of Fluid Mechanics, vol. 807, 2016, pp. 155-166.
[29]. Е.М. Миронова. Математическое моделирование движения водных потоков, снежных лавин и селей. Диссертация на соискание ученой степени кандидата физико-математических наук. МГУ, 1987/ E. Mironova. Mathematical Modelling of the Motion of Water Flows, Snow Avalanches, and Floods. PhD thesis, Lomonosov Moscow State University, 1987 (in Russian).
[30]. Kenichi Oda, Shuji Moriguchi et al. Simulation of a snow avalanche model test using computational fluid dynamics. Annals of Glaciology, vol. 52, issue 58, 2011, pp. 57-64.
[31]. Д.И. Романова Трёхмерное моделирование потоков жидкости Хершеля-Балкли на склоне в OpenFOAM. Труды ИСП РАН, том 29, вып. 1, 2017 г., стр. 85-100 / Romanova D.I. 3D flow modeling of Herschel-Bulkley fluid on the slope in OpenFOAM. Trudy ISP RAN/Proc. ISP RAS, vol. 29, issue 1, 2017, pp. 85-100 (in Russian). DOI: 10.15514/ISPRAS-2017-29(1)-6.
[32]. Daria Romanova. Comparison of single-velocity and multi-velocity multiphase models for slope flow simulations. In Proc. of the 2020 Ivannikov Ispras Open Conference (ISPRAS), 2020, pp. 170-174.
[33]. Д.И. Романова. Архитектура программного средства с открытым исходным кодом для численного моделирования потоков на горных склонах. Труды ИСП РАН, том 32, вып. 6, 2020 г., стр. 183-200
/ D.I. Romanova. Design of open source software architecture for numerical modeling of flows on mountain slopes. Trudy ISP RAN/Proc. ISP RAS, vol. 32, issue 6, 2020, pp. 183-200 (in Russian). DOI: 10.15514/ISPRAS-2020-32(6)-14.
[34]. P. Sampl. Current status of the avl avalanche simulation model—numerical simulation of dry snow avalanches. In Proc. of the "Pierre Beghin" International Workshop on Rapid Gravitational Mass Movements, 1993, pp. 269-296.
[35]. S.B. Savage, K. Hutter. The motion of a finite mass of granular material down a rough incline. Journal of Fluid Mechanics, vol. 199, 1989, pp. 177-215.
[36]. S. B. Savage, K. Hutter. The dynamics of avalanches of granular materials from initiation to runout. part i: Analysis. Acta Mechanica, vol. 86, issue 1, 1991, pp. 201-223.
[37]. Thomas Scheiwiller. Dynamics of powder-snow avalanches. PhD thesis, ETH Zurich, 1986.
[38]. Sherif H. El Tahry. k-epsilon equation for compressible reciprocating engine flows. Journal of Energy, vol. 7, issue 4, 1983, pp. 345-353.
[39]. A. Voellmy, A. Roch. Über die Zerstörungskraft von Lawinen. Schweizerische Bauzeitung -Wochenschrift für Architektur, Ingenieurwesen und Maschinentechnik, 73, 1955 (in German).
[40]. Н.А. Володичева, Е.М. Миронова др. Использование математического моделирования для определения границ распространения лавин. Материалы гляциологических исследований, вып. 56, 1986 г., стр. 78-81 / N. Volodicheva, E. Mironova et al. The use of mathematical modelling to determine the boundaries of propagation of avalanches. Materialy Glyatsiologicheskikh Issledovaniy, issue 56, 1986, pp. 78-81 (in Russian).
[41]. Yuya Yamaguchi, Shinsuke Takase et al. Three-dimensional nonstructural finite element analysis of snow avalanche using non-newtonian fluid model. Transactions of the Japan Society for Computational Engineering and Science, 2017, Paper No. 20170011 (in Japanese).
Информация об авторах / Information about authors
Дарья Игоревна РОМАНОВА получила степень магистра на кафедре гидромеханики механико-математического факультета МГУ в 2017 году, в настоящее является младшим научным сотрудником лаборатории вычислительных методов механико-математического факультета МГУ и стажёром-исследователем ИСП РАН. Она является разработчиком программного обеспечения с открытым исходным кодом для численного моделирования задач механики сплошных сред, включая турбулентные течения, течения со свободной поверхностью, многофазные течения с фазовым переходом и включениями твердых частиц.
Daria Igorevna ROMANOVA received her master's degree at the Department of Hydromechanics, Faculty of Mechanics and Mathematics, Moscow State University in 2017, and is currently a junior researcher at the Laboratory of Computational Methods at the Faculty of Mechanics and Mathematics of Moscow State University and an intern-researcher at ISP RAS. She is an open source software developer for the numerical modeling of continuum mechanics problems, including turbulent flows, free surface flows, multiphase flows with a phase transition, and solid inclusions.
