Трёхмерное моделирование схода лавинных потоков средствами пакета OpenFOAM Текст научной статьи по специальности «Физика»

Трёхмерное моделирование схода лавинных потоков средствами пакета

OpenFOAM

Д.И. Романова <romanovadi@gmail.com> Московский государственный университет имени М.В. Ломоносова, 119991, Россия, Москва, Ленинские горы, д. 1 Федеральное государственное учреждение "Федеральный научный центр Научно-исследовательский институт системных исследований Российской

академии наук", 117218, Москва, Нахимовский просп., 36, к.1

Аннотация. В работе создана модель схода лавинного потока в пакете OpenFOAM, где лавина смоделирована как двухфазный турбулентный поток (снег - воздух). За основу был взят численный метод для решения задач движения двухфазных сред с поверхностью раздела - метод переноса объёмной доли (VOF). Данный алгоритм реализован в решателе InterFoam, который и был использован в работе. Использовалась К - s модель турбулентности. В модели снег был представлен как нелинейно вязкая жидкость, описываемая реологическими соотношениями Хершеля-Балкли. Воздух представляет собой вязкую ньютоновскую среду. Для описания модели используются осреднённые по Рейнольдсу уравнения Навье-Стокса, реологические соотношения, а также уравнения для турбулентной кинетической энергии и диссипации. Построена расчётная область лавинного очага номер 22 горы Юкспор Хибинских гор на основе цифровой модели рельефа формата ASCII GRID. Взяты карты зоны зарождения и зоны отложений 129-ой лавины, глубина снежного покрова в зоне зарождения составила 1,5 метра, рассчитан поток, возникающий в зоне зарождения, и сравнивается модельная зона отложения с натурными данными. Также смоделировано движение снего-пылевого облака. Получены распределения скоростей, давления и объёмной доли снега во все моменты времени, форма лавинных отложений. Средняя скорость потока составила 44,8 м/с, что близко к действительным скоростям движения лавин в данном очаге. Максимальная скорость лавинного потока (включая снего-пылевое облако) составила 78 м/с. Данные результаты позволяют оптимально проектировать противолавинные сооружения и определять лавиноопасные территории.

Ключевые слова: лавиноведение; турбулентный поток; двухфазная среда; жидкость Хершеля-Балкли.

DOI: 10.15 514Л SPRAS-2016-29( 1 )-6

Для цитирования: Романова Д.И. Трёхмерное моделирование потоков жидкости Хершеля-Балкли на склоне в OpenFOAM. Труды ИСП РАН, том 29, вып. 1, 2017 г., стр. 85-100. DOI: 10.1551MSPRAS-2017-29(l)-6

1. Введение

Склоновые потоки, в том числе лавины, — распространённое явление в горах. В России они встречаются на Кавказе, в Хибинах, горах Прибайкалья, Забайкалья, на Урале, Чукотке — словом, всюду, где углы падения склонов больше 15°, и имеется снежный покров глубиной 30 - 40 см и выше. Громадные разрушения, вызываемые обвалами снега, человеческие жертвы, миллионные убытки, невозможность поддержать круглогодичное бесперебойное движение по горным дорогам, требуют уделять большое внимание проблеме борьбы с лавинами, исследовать и моделировать их. Для оценки лавинной опасности при изысканиях и строительстве в лавиноопасных районах необходимо определять наибольшую дальность выброса лавин и возможную силу удара лавинного снега. Для определения силы удара требуется знать скорость движения лавины (точнее, её фронтальной части), а также высоту фронта.

2. Модель движущейся среды

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

Уравнения Рейнольдса выглядят следующим образом (черта над буквой означает осреднение по Рейнольдсу):

д(рщ) д (_-- —;-гл др дтij

— + ~j (диЪ + Pu >J) =P9i-—j + TxJ>

где f у - компоненты осреднённого по Рейнольдсу тензора вязких напряжений:

fy = 2/хёу, ц = n(\f\),

+ 1Я=л/2Ш- h(e) = eijei.

Далее ц выражается через кинематическую вязкость: /./ = pv, а р и v вычисляются через объёмную долю снега а v = (1 - a)vair + avsnow, р = (1 - a)pair + apsnow,

!1 — снег 0 < а < 1 — в переходном слое 0 — воздух

Здесь индекс «air» соответствует воздуху, а индекс «snow» - снегу.

Для различных фаз использованы различные реологические соотношения, для снега взяты реологические соотношения Хершеля-Балкли [1, 2], воздух представлен как ньютоновская среда.

Ньютоновская модель: Модель Хершеля-Балкли:

т0+к\у\1<

т = цу, ц = const Здесь т - сдвиговое напряжение.

т = цу, ¡1

m

к,п = const

■ скорость сдвига.

В случае простого сдвигового потока соответствующие реологические соотношения имеют вид (см. рис. 1)

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

Fig. 1. The shear rate dependency of shear stress in a simple shear flow for different

rheological models.

В силу выбранных реологических соотношений кинематическая вязкость вычисляется по-разному для воздуха и для снега:

5 „2/--------,, TQ+fcl/Г

vair = const = 1.48 ■ 10

м2/сек,

v

= min ^v.

SUOWQ>

Psnowl

где vsnow = const и задаётся при калибровке модели в соответствии с данными исследуемой лавины. Необходимость введения константы v.

SUOWq

связана с тем, что значение

кинематической вязкости снега при т0 Ф 0, стремится к бесконечности, когда \у\ -> 0, как показано на рис. 2.

В расчетах задана плотность воздуха pair = 1 кг/м3. Параметры для снега, такие как psnow>k,n,T0, будут подобраны при калибровке модели в соответствии с натурными данными. Скорости фаз считаются одинаковыми.

Y

Puc. 2. Зависимость эффективной вязкости от скорости сдвига для среды Хершеля-

Бапкпи при п < 1.

Fig. 2. The shear rate dependency of effective viscosity' at n < 1.

Для описания модели использованы следующие уравнения. Уравнение неразрывности:

£Р _l_ d(puj) _ Q dt dxt

Уравнение объёмной доли фазы:

За d(auj) _ ^ dt dxj ~

В качестве модели турбулентности берётся К — £ модель. Уравнение для турбулентной кинетической энергии К:

д(Рк) d(pujK) _д(ак\ а (ßtdK\

dt + dxj ~ dxj^dxj) + dxj \aKdxj) P£-

Уравнение для диссипации e:

l=c p i_ с Ei + J-itLlLl

dt dxj E1 к К И EZ К dx]\aEdx]J'

-pu\u'j = 2ntäi} -\р8цК, nt = pCp ^ PK = -pu\u'j Щ

Приведённая вьше модель турбулентности содержит пять констант [1], которые задаются следующим образом:

= 0.09; СЕ1 = 1.44; Се2 = 1.92; ок = 1.0; оЕ = 1.3.

3. Программное обеспечение

Модель лавинного потока была построена в пакете OpenFOAM. OpenFOAM (Open Source Field Operation And Manipulation CFD ToolBox) -открытая интегрируемая платформа для численного моделирования задач механики сплошной среды. Используется решатель InterFoam, предназначенный для расчета нестационарного течения двух сред, разделённых 88

границей раздела или свободной поверхностью. Для нахождения формы свободной поверхности двухфазный алгоритм в ПисгРоат использует один из методов, в котором интерфейс не рассматривается как резкая граница (метод захвата интерфейса). Расчет происходит на фиксированной сетке, включающей свободную поверхность. Форма свободной поверхности определяется расчётом доли заполненности каждой ячейки вблизи поверхности. Это можно сделать, решая уравнение транспорта для доли ячейки, занимаемой жидкой фазой; такой метод называется "объём жидкости"(УоИшге-оНЧшс!) или схема УОБ. Более подробное описание использованного решателя можно найти в [3].

4. Объект моделирования

Схема лавинных очагов в районе г. Кировека

Рис. 3. Схема лавинных очагов в районе г. Кировека. Fig. 3. Scheme of avalanche sites nearby city Kirovsk

При помощи построенной модели была рассчитана лавина, сошедшая в 1965 году из 22-ого лавинного очага горы Юкспор в Хибинах, вблизи города Кировск, см. рис. 3,4.

Именно в данном очаге был произведён профилактический спуск лавины в 22:00 18 февраля 2016 года. Лавина приобрела катастрофический характер и унесла жизни троих людей, в том числе одного из сотрудников лавинной службы. Были засыпаны железная и автомобильная дороги, выбиты стёкла в трёх близлежащих домах.

На основе цифровой растровой карты 22-ого лавинного очага формата ASCII GRID была построена расчётная область (рис. 5), которая представляет собой пространственную область, ограниченную снизу поверхностью склона, а сверху- поверхностью такой же формы, сдвинутой вверх на 25 метров. Из паспорта лавины, сошедшей в 1965г., были взяты следующие параметры:

• максимальная толщина слоя снега в месте отрыва лавины -1.5 м;

• средняя толщина слоя снега в месте отрыва лавины - 0.8 м;

• максимальная плотность лавинных отложений - 320 кг/мЗ;

• средняя плотность лавинных отложений - 290 кг/мЗ;

• минимальная плотность лавинных отложений - 260 кг/мЗ;

• состояние снега - сухой;

• максимальная толщина лавинных отложений -4 м;

Рис. 4. Фотография лавинного очага №22 (НИЛ снежных лавин и селей)

Fig. 4. Photo of 22 avalanche cite (Research Laboratory of Snow Avalanches and

Debris Flows).

• средняя толщина лавинных отложений - 0.72 м.

X

Рис. 5. Расчётная область. Fig. 5. Computational domain.

5. Расчётная сетка

Расчёт производился на сетке из ячеек в форме параллелепипедов с линейным размером 5м, состоящей из:

• 206 630 точек

• 532 887 граней

• 447 345 внутренних граней

• 163 372 ячеек

Puc. 6. Сеточные линии расчётной области. Fig. б. Mesh grid of computational domain.

6. Начальные и граничные данные

В задаче заданы следующие граничные условия.

На нижней границе расчётной области установлено условие прилипания. Боковые и верхняя границы расчётной области заданы как стандартные участки границы с условием протекания, на которых действует атмосферное давление. На рис. 7 зелёным цветом обозначен слой снега, лежащий в зоне зарождения лавины, при t = 0 он покоится (U = 0 всюду). Глубина снега во всех точках при t = 0 не превышает 1.5 м. Математически это выражается так, что маркерная функция а равна 1 в тех ячейках расчётной сетки, где снег есть, и 0 в остальных (а = a(x,y,z)). Также при t = 0: г = 0, К = 0, ßt = 0.

Серым цветом на рис. 7 изображена реальная зона отложений исследуемой лавины.

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

X

Рис. 7. Лавинный поток в начальный момент времени. Fig. 7. Avalanche at initial time.

7. Расчёт

В результате дискретизации исходных уравнений в частных производных методом конечных объёмов были получены системы линейных алгебраических уравнений, которые решались итерационными методами.

Для решения дискретизированных уравнений для объёмной доли фазы а, турбулентной кинетической энергии К, диссипации е и уравнения движения используется решатель со сглаживанием с помощью метода Гаусса-Зейделя. Уравнение на давление, которое получено при помощи уравнения неразрывности, решается методом сопряжённых градиентов с предобуславливателем.

Расчёты производились на высокопроизводительном вычислительном сервере НИИСИ РАН, состоящем из:

• 2 х Xeon E5-2670vl (Sandy Bridge) 8 cores 2.6 GHz (16 threads)

• 4xchannel DDR3 memory controller

• MEM 256 GB DDR3-1333 ECC

• OS Ubuntu 15.10 64-bit

8. Результаты

В процессе вычислений были получены значения щ(х,у,z, t) i = 1,2,3, a(x,y,z,t), p(x,y,z,t), а также граница лавинного потока, которая была определена из условия a = 0.1.

Проведены расчёты со значениями параметров psnow, vvn0W(| , т0, к, п для снега в следующих диапазонах:

р [кг/м3] v0 [м2/с] Т°/р К/с2] к/р [м2/с] п

200 100 1 10"5 2

300 105 0.1 10"4 0.5

10"2 10 10"6

Было посчитано 108 случаев, что заняло 36 часов машинного времени. Были сделаны выводы, что параметры т0 и v0 влияют на динамику потока, и при подходящем выборе поток начинает замедляться на склоне. Скорости потока в данных калибровочных расчётах получились в среднем 60 м/с, что много для данной лавины. Было найдено, что на величину скорости существенно влияет параметр вязкости к и проведено дополнительное исследование по параметру к с остальными зафиксированными подходящими параметрами, взятыми из предыдущих расчётов. Рассчитывались варианты с к/р = 2, 3, 4, 5, 6, 7, 8, 9, 10, 50, 100 м2/с, Т°/р = Ю м2/с2, v0 = 105 м2/с, р = 200 кг/м3, их — 0.5. Расчёт с параметром ^/р = 5м2/с наиболее точно описал натурные данные.

На рис. 8-11 изображены границы лавинного потока и распределение величины скорости на поверхности потока в разные моменты времени. На рис. 12 - 15 изображено распределение величины скорости на верхней границе расчётной области в те же моменты времени. Оно отображает распространение снего-пылевого облака.

Оси х, у - горизонтальные, ось z - вертикальная. На цветной легенде U Magnitude - модуль скорости, U Magnitude = ^и2 + Uy + .

Time: 0.000000

U Magnitude

_1.530е+02

Рис. 8. Лавинный поток в начальный момент времени. Fig. 8. Avalanche at start moment.

U Magnitude

г 1,530e+02

Time: 25.000000

Рис. 9. Ускорение лавинного потока. Fig. 9. Avalanche acceleration

U Magnitude

,-1.530e402

!

Time: 35.000000

u

U Magnitude

_-1.55Ce+C2

Puc. 10. Замедление лавинного потока. Fig. 10. Avalanche deceleration.

Time: 45.000000

Puc. 11. Остановка лавинного потока. Fig. 11. Avalanche stop.

U Magnitude

r 1.417В+02

Time: 0.000000

U Magnitude

Рис. 12. Снего-пыпевое облако в начальный момент времени. Fig. 12. Snow-dust cloud at start moment.

Time: 25.000000

Рис. 13. Ускорение снего-пыпевого облака. Fig. 13. Snow-dust cloud acceleration

U Magnitude

_ö.l50e+01

Time: 35.000000

Рис. 15. Замедление сне го-пылевого облака. Fig. 15. Snow-dust cloud deceleration.

U Magnitude

_-й.<?1йе4Ш

Рис. 14. Движение снего-пыпевого облака. Fig. 14. Snow-dust cloud motion.

Time: 45.000000

Можно видеть, что в данном расчёте после полной остановки потока форма лавинных отложений, полученная в процессе вычислений, близка к натурной форме лавинных отложений. Средняя скорость потока составляет 44,8 м/с, что близко к действительным скоростям движения лавин в данном очаге. Максимальная скорость лавинного потока (включая снего-пылевое облако) составляет 78 м/с.

9. Заключение

Была создана модель лавинного потока в пакете OpenFOAM, позволяющая получить значения компонент скорости, давление, объёмную долю снега в каждой точке расчётной области в различные моменты времени. Была рассчитана лавина, сошедшая в Хибинах, на горе Юкспор, в 22-ом лавинном очаге. По полученным данным была построена зона лавинных отложений, границы которой близки к зафиксированным непосредственно после схода реальной лавины. Вычислены распределения скоростей, объёмной доли снега и давления во все моменты времени, в частности максимальная скорость лавинного потока. Данные результаты позволяют определить лавиноопасную зону и помогают оптимально спроектировать защитные сооружения. В дальнейшем планируется рассчитать задачу на более подробной сетке, а также построить для сравнения данную модель в другом пакете, например, INMOST или при помощи собственного кода.

Автор благодарит профессора МГУ М. Э. Эглит за многочисленные обсуждения модели, П.Б. Богданова, сотрудника ФГУ ФНЦ НИИСИ РАН, за помощь в организации вычислительного алгоритма и его программной реализации и Ю.Г. Селиверстова, Т.Г. Глазовскую, A.C. Турчанинову сотрудников Научно-исследовательской лаборатории снежных лавин и селей Географического факультета МГУ за предоставленные данные и обсуждение результатов работы.

Список литературы

[1]. Эглит М.Э., Якубеико А.Е. Влияние захвата донного материала и неньютоновской реологии на динамику турбулентных склоновых потоков. Известия Российской академии наук. Механика жидкости и газа. 2016, № 3, стр. 3-15. DOI: 10.7868/S056852811603004X

[2]. Eglit М.Е., Yakubenko А.Е. Numerical modeling of slope flows entraining bottom material. Cold Regions Science and Technology 108 (2014), 139-148

[3]. Joel H. Ferziger, Milovan Perit. Computational Methods for Fluid Dynamics - 3., rev. ed. - Berlin; Heidelberg; New York; Barcelona; Hong Kong; London; Milan; Paris; Tokyo: Springer, 2002.

3D avalanche flow modeling using OpenFOAM

D.I. Romanova<romanovadi@gmail. com> Lomonosov Moscow State University, GSP-1, Leninskie Gory, Moscow, 119991, Russia. Federal State Institution «Scientific Research Institute for System Analysis of the Russian Academy of

Sciences»,

36k1, Nakhimovskiy pr., Moscow, 117218, Russia

Abstract. In this paper, a model of a snow avalanche was created with the help of open source CFD software OpenFOAM. The avalanche is considered as a turbulent two-phase flow — snow and air. We take incompressible Herschel-Balkley fluid as a model of snow. Air is a newtonian fluid. For tracking and locating the free surface we use the volume of fluid (VOF) method. In calculations we use solver interFoam which is based on the VOF method. The K - s turbulence model was used. Navier-Stokes equations, rheological ratios, equations for turbulent kinetic energy and dissipation are used to determine the model. The avalanche occurred at the 22nd site on the Ukspor mountain was modeled. Computational domain was made using digital terrain model in ASCII GRID format. The shape of snow deposits area was calculated and compared with the real data. The velocity field of flow, pressure distribution, the field of volume snow fraction were obtained for different time instances. The value of the average velocity of the flow was 44,8 m/s, the value of the maximum velocity of the flow (including snow-dust cloud) was 78 m/s. These results allow to determine the avalanche hazard area and help to optimally design the defense systems. In future it is planned to simulate this event using more detailed grids and realize this model in other software for example INMOST or with own code to compare the results.

Keywords: glaciology; turbulent flow; two-phase medium; Herschel-Bulkley fluid. DOI: 10.15514/ISPRAS-2017-29(l)-6

For citation: 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. DOI: 10.15514/ISPRAS-2017-29(1 )-6

References

[ 1 ]. Eglit M.E., Yakubenko A.E. Effect of the bottom material capture and the non-Newtonian rheology on the dynamics of turbulent downslope flows. Fluid Dynamics. 2016, vol. 51, №3, pp. 299-310. DOI: 10.1134/S0015462816030017

[2]. Eglit M.E., Yakubenko A.E. Numerical modeling of slope flows entraining bottom material. Cold Regions Science and Technology, 108 (2014), pp. 139-148

[3]. Joel H. Ferziger, Milovan Perit. Computational Methods for Fluid Dynamics - 3., rev. ed. - Berlin; Heidelberg; New York; Barcelona; Hong Kong; London; Milan; Paris; Tokyo: Springer, 2002