CC BY
33
УДК 551.4.013 DOI 10.18522/0321-3005-2016-3-20-24
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ДВИЖЕНИЯ ОБВАЛОВ С ИСПОЛЬЗОВАНИЕМ КОНТИНУАЛЬНОГО ПОДХОДА*
© 2016 г. Н.С. Орлова, М.В. Волик
Орлова Наталья Сергеевна - кандидат технических наук, старший преподаватель, кафедра математики и информатики, Финансовый университет при Правительстве РФ; ученый секретарь, Южный математический институт - филиал Владикавказского научного центра Российской академии наук, ул. Маркуса, 22, г. Владикавказ, 362027, e-mail: norlova.umi.vnc@gmail.com
Волик Мария Владимировна - старший преподаватель, кафедра математики и информатики, Финансовый университет при Правительстве РФ; младший научный сотрудник, Южный математический институт - филиал Владикавказского научного центра Российской академии наук, ул. Маркуса, 22, г. Владикавказ, 362027, e-mail: volikmv@mail.ru
Orlova Natalya Sergeevna - Candidate of Technical Science, Senior Lecturer, Department of Higher Mathematics, Financial University under the Government of RF; Scientific Secretary, Southern Mathematical Institute of VSC RAS, Markus St., 22, Vladikavkaz, 362027, Russia, e-mail: norlova.umi.vnc@gmail.com
Volik Maria Vladimirovna - Senior Lecturer, Department of Higher Mathematics, Financial University under the Government of RF; Junior Researcher, Southern Mathematical Institute of VSC RAS, Markus St., 22, Vladikavkaz, 362027, Russia, e-mail: volikmv@mail.ru
Рассматривается движение обвала по склону, сопряженному с горизонтальным участком. Теоретическое исследование осуществлялось с использованием свободно распространяемого пакета для решения прикладных задач гидро- и аэромеханики OpenFOAM. Использовалась двухжидкостная модель на основе континуального подхода и кинетической теории гранулярного газа. Представлены результаты двумерных численных расчетов распределения обвальной массы после ее падения в зависимости от угла склона к горизонтальной поверхности.
Ключевые слова: математическое моделирование, обвал, ожижение, континуальный подход, кинетическая теория гранулярных газов, пакет OpenFOAM, решатель twoPhaseEulerFoam.
The rockfall movement of the slope, which is associated with the horizontal section, was investigated. The theoretical investigation was performed using redistributable package OpenFOAMfor solution of applied problems of hydro and aerodynamics. Continuum approach and kinetic theory of granular gas were used. Calculations were performed for various values of the slope angle to horizontal surface. Two-dimensional numerical calculations of distribution of rockfall mass were presented.
Keywords: mathematical modeling, rockfall, fluidization, continuum approach, kinetic theory of granular gas, OpenFOAM, twoPhaseEulerFoam.
В настоящее время оценка размеров зон поражения, вызванных обвалами массы горных пород, представляется актуальным научно-практическим исследованием. Осуществить его можно с использованием математического моделирования. Несмотря на достаточно большое количество известных зарубежных и отечественных научных работ в этом направлении, отсутствует математическая модель обвалов, которая бы учитывала ожижение обломков, вызванное их хаотическим движением. Под ожижением понимается изменение объемной плотности обломков по толщине слоя, которая может меняться во время его движения.
Для описания обвалов используются в основном дискретные модели, которые описывают движение потока вещества в виде движения совокупности отдельных структурных частиц. Для каждой вы-
бранной структурной частицы вещества определяются ее динамические характеристики (текущие координаты, вектор скорости) на каждом последующем временном шаге с учетом ее параметров и параметров контактирующих с ней частиц на предыдущем шаге. Дискретные модели могут описывать движение потока вещества либо в виде группы отдельных не соударяющихся друг с другом частиц, каждая из которых без сопротивления проходит сквозь другую (соседнюю) частицу, либо в виде группы упруго соударяющихся друг с другом частиц [1—5]. Таким образом, дискретные модели позволяют моделировать ограниченное количество движущихся обломков, увеличение количества которых требует достаточно мощных вычислительных ресурсов, так как необходимо моделировать движение каждого обломка в отдельности.
*Работа выполнена при поддержке гранта РФФИ № 16-35-00147.
Известны и непрерывные модели обвалов [6, 7], в которых движение вещества представляется в виде сплошной среды, характеризующейся неразрывным полем значений физических параметров (скорости, давления, сил). Но эти модели не учитывают процесс ожижения потока обвальных горных пород, движущихся по склону.
В данной работе исследуется математическая модель движения обломков горных пород с использованием континуального подхода, когда движение потока вещества представляется в виде сплошной среды. Такой подход не требует использования мощных вычислительных ресурсов. Кроме того, в модели учитывается ожижение потока обвальных горных пород в процессе их движения по склону. Моделирование ожижения потока обвальных горных пород осуществляется с использованием кинетической теории гранулярных газов, которая учитывает хаотическое движение обломков (вследствие их столкновений друг с другом и с поверхностью склона) как в плотном, так и в разреженном состоянии.
Таким образом, кинетическая теория гранулярных газов дает возможность учесть влияние столкновений обломков друг с другом и со склоном на движение обвалов. Ранее такой подход не применялся для описания обвалов.
В данной работе для моделирования движения обвалов использовался свободно распространяемый пакет для численного моделирования задач гидроаэромеханики OpenFOAM (англ. Open Source Field Operation And Manipulation CFD ToolBox) при поддержке программы «Университетский кластер» с удаленным доступом к консоли на управляющем узле вычислительного кластера JSCC RAS веб-лаборатории UniHUB (URL: https://unihub.ru/ tools/js3console (date of treatment: 30.04.2016)).
Для описания движения обвала с учетом ожижения использовался решатель twoPhaseEulerFoam, в котором реализована двухжидкостная модель кипящего (ожиженного) слоя на основе континуального подхода (подхода Эйлера), при котором движение слоя рассматривается как движение двух взаимодействующих континуумов, связанных с газом и частицами. Основные уравнения двухжидкостной модели - уравнения неразрывности и уравнения количества движения для обеих фаз [8-11].
-|(афРф)+^-(афрфиф)= 0; (1)
-|(a9P9U9)+V-(a9p9U9U9)+V-(a9R )=
= -афУр + aфPфg + М ф. (2)
Индекс ф означает принадлежность к фазе (твердой a или газовой b, обе фазы считаются несжимаемыми); аф - объемная доля соответствующей фазы; Рф, Пф - плотность и вектор скорости фазы; R^eff - тензор эффективных напряжений; P -давление газовой фазы; g - ускорение свободного падения; Мф - член, моделирующий обмен импульсом между фазами. Выражения для коэффициентов и членов, входящих в уравнения (1), (2), подробно описаны в литературе [8-10].
Для учета эффектов, обусловленных взаимодействием частиц друг с другом, используется кинетическая теория (по аналогии с кинетической теорией газа), с помощью которой можно выразить эффективные напряжения, возникающие в дисперсной фазе за счет движения и столкновений частиц друг с другом. По аналогии с термодинамической температурой вводится гранулярная температура 0 как средняя энергия флуктуаций скорости частиц (обломков). В работе [9] приведено уравнение для расчета гранулярной температуры.
Проведены тестовые расчеты модели с целью дальнейшего ее использования для моделирования движения обвалов.
Задача решалась в двумерном приближении. Использовались следующие начальные условия:
Па = 0,25 ; Пь = 0; P = 0; аа = 0,6 ; 0 = 0 .
На рис. 1 представлена область расчетов (рис. 1а) и граничные условия задачи (рис. 1б). Рассматривалась область, состоящая из склона, сопряженного с горизонтальным участком. Расчеты проводились при разном значении угла склона к горизонтальной поверхности (Pi = 37°, р2 = 45° и р3 = 54°).
Размеры вычислительной области: высота -180 м; ширина - 140 м. Движение обвала рассчитывалось за 1300 с. Это позволило рассмотреть движение обвальной массы по склону и горизонтальной поверхности до момента, когда она остановится. При этом использовался шаг по времени, равный 0,005 с. Шаг по пространству равен 1 м. Проводилось распараллеливание расчетов на 8 ядрах.
Далее подробно представлены граничные условия на стенках (wall):
dP . 9а а п 90 „
Па = 0 ; иь = 0; — = 0;—а = 0 ; — = 0 ,
а ' ь 'л ''i 'л '
9n on 9n а также на входной левой границе (inlet):
9P 90
Па = 0 ; Пь = 0; — = 0; аа = 0; — = 0 .
а ? ь 'а 'л
on 9n
Граничные условия на верхней и правой свободных границах (outlet):
дЦ^ = 0 .dU,
дР
да„
дв
дп
— = 0 ; — = 0;даа = 0 ; — = 0 . дп дп дп дп
Рис. 1. Область расчетов
В верхней части склона в начальный момент времени располагалась обвальная масса в виде треугольной области площадью сечения примерно 38 м2, как это видно на рис. 1б. Обвальной массе в начальный момент времени задавалась скорость 0,25 м/с. Обвальная масса представляется как совокупность твердых частиц сферической формы. Диаметр частиц (обломков) принимался одинаковым и равным 5 мм. В дальнейшем планируется увеличивать размер частиц (обломков). Следует отметить, что увеличение размера частиц предполагает увеличение пространственного шага в рамках континуального подхода.
На рис. 2 представлены графики распределения объемной доли обвальной массы (рис. 2а) и ее скорости иа (рис. 2б) по нормали к наклонной поверх-
ности в точке с координатами х = 18 м, у = 16 м в момент времени t = 300 с при в1 = 37°. Результаты расчетов на рис. 2 показывают неравномерное распределение обвальной массы, а также ее скорости по высоте слоя. Это свидетельствует об ожижении обвальной массы в процессе ее движения по склону.
На рис. 3 отражена оценка распределения толщины слоя Н вдоль горизонтального участка после падения обвальной массы ^ = 1000 с). Кривая 1 на рис. 3 соответствует в = р1; 2 - в = р2; 3 - в = р3.
Из рис. 3 видно, что с увеличением угла склона к горизонтальной поверхности в максимальное значение толщины слоя увеличивается, а расстояние, пройденное обвальной массой вдоль горизонтальной поверхности, уменьшается (при в1 ~ 55 м; при в2 ~ 50; при в3 ~ 40 м).
alpha
0,20
0,15
0,10
0 0,5 1,0 1,5 2,0 2,5 3,0 3,5 4,0 4,5 5,0 П - ,V1
Оа, М/С
0,10
0,04
0,02
VGI
0 0,50 1,0 1,5 2,0 2,5 3,0 3,5 4,0 4,5 5,0
п, м
VGI
Рис. 2. Распределение обвальной массы и ее скорости по нормали к наклонной поверхности
Рис. 3. Распределение толщины слоя вдоль горизонтального участка
0,30
0,25
0,05
0
0,14
0,12
0,0«
0,06
0
Проведено сравнение полученных расчетов с результатами для других известных моделей, которые описывают реальные наблюдения. Полученные размеры зоны поражения в среднем совпадают с результатами, которые рассчитаны по упрощенным динамическим моделям и при обработке данных наблюдений [3, 12].
Таким образом, можно сделать вывод о том, что данную модель можно использовать для описания движения обвалов и оценки зон поражения при обвалах.
Литература
1. Михайлов В.О. Классификация численных математиче-
ских моделей селевых и склоновых процессов // Инженерная геология. 2011. № 3. С. 56-63.
2. Михайлов В.О. Трехмерная математическая модель об-
вальных процессов // Вестн. Моск. ун-та. Сер. 5: География. 2011. № 4. С. 53-58.
3. Luuk K.A. Dorren A review of rockfall mechanics and mod-
eling approaches // Progress in Physical Geography. 2003. Vol. 27, № 1. P. 69-87.
4. Rammer W., Brauner M., Dorren L.K.A., Berger F.,
Lexer M.J. Evaluation of a 3-D rockfall module within a forest patch model // Natural Hazards ahd Earth System Sciences. 2010. № 10. P. 669-711.
5. Hungr O. A model for the runout analysis of rapid flow
slides, debris flows and avalanches // Canadian Geotech-nical J. 1995. Vol. 32, № 4. P. 610-623.
6. Божинский А.Н., Назаров А.Н. Динамика двухфазного
селевого потока // Вестн. Моск. ун-та. Сер. 5: География. 1999. № 5. С. 15-20.
7. Fangqiang Wei, Hu Kaiheng, Jose Luis Lopez, Cui Peng.
Method and its application of the momentum model for debris flow risk zoning // Chinese Science Bulletin. 2003. № 48(6). Р. 594-598.
8. Rusche H. Computational Fluid Dynamics of Dispersed
Two-Phase Flows at High Phase Fractions: Thesis submitted for the degree of Doctor of Philosophy of the University of London and Diploma of Imperial College. L., 2002. 343 p.
9. Wachem B. van. Derivation, Implementation, and Validation
of Computer Simulation Models for Gas-Solid Fluidized Beds: Dissertation at Delft University of Technology. Delft, 2000. 222 p.
10. Каменецкий Е.С., Орлова Н.С., ВоликМ.В., Минасян Д.Г.
Исследование динамики кипящего гранулированного слоя с использованием пакета OpenFOAM // Изв. вузов. Сев.-Кавк. регион. Естеств. науки. 2014. № 5. С. 37-42.
11. Gidaspow D. Multiphase flow and fluidization: Continuum
and kinetic theory descriptions. Boston, 1994. 211 p.
12. Nilsen M. Modelling of rockfall runout range. Employing
empirical and dynamical methods: Master Thesis in Geo-
Поступила в редакцию
sciences Discipline: Environmental geology and geo-hazards. Oslo, 2008. 96 p.
References
1. Mikhailov V.O. Klassifikatsiya chislennykh mate-
maticheskikh modelei selevykh i sklonovykh protsessov [Classification of numerical mathematical models of mud-flow and slope processes]. Inzhenernaya geologiya, 2011, no 3, pp. 56-63.
2. Mikhailov V.O. Trekhmernaya matematicheskaya model'
obval'nykh protsessov [A three-dimensional mathematical model of the landslide processes]. Vestn. Mosk. un-ta. Seriya 5. Geografiya, 2011, no 4, pp. 53-58.
3. Luuk K.A. Dorren A review of rockfall mechanics and
modeling approaches. Progress in Physical Geography, 2003, vol. 27, no 1, pp. 69-87.
4. Rammer W., Brauner M., Dorren L.K.A., Berger F., Lexer
M.J. Evaluation of a 3-D rockfall module within a forest patch model. Natural Hazards and Earth system Sciences, 2010, no 10, pp. 669-711.
5. Hungr O. A model for the runout analysis of rapid flow
slides, debris flows and avalanches. Canadian Geotech-nical J, 1995, vol. 32, no 4, pp. 610-623.
6. Bozhinskii A.N., Nazarov A.N. Dinamika dvukhfaznogo
selevogo potoka [The dynamics of a two-phase debris flow]. Vestn. Mosk. un-ta. Seriya 5. Geografiya, 1999, no 5, pp. 15-20.
7. Fangqiang Wei, Hu Kaiheng, Jose Luis Lopez, Cui Peng.
Method and its application of the momentum model for debris flow risk zoning. Chinese Science Bulletin, 2003, no 48(6), pp. 594-598.
8. Rusche H. Computational fluid dynamics of dispersed two-phase flows at high phase fractions. Thesis submitted for the degree of Doctor of Philosophy of the University of London and Diploma of Imperial College. London, 2002, 343 p.
9. Wachem B. van. Derivation, implementation, and validation
of computer simulation models for gas-solid fluidized beds. Dissertation at Delft University of Technology. Delft, 2000, 222 p.
10. Kamenetskii E.S., Orlova N.S., Volik M.V., Minasyan D.G.
Issledovanie dinamiki kipyashchego granulirovannogo sloya s ispol'zovaniem paketa OpenFOAM. [The study of the dynamics of fluidized granular layer using the package OpenFoam]. Izv. vuzov. Sev.-Kavk. region. Estestv. nauki, 2014, no 5, pp. 37-42.
11. Gidaspow D. Multiphase flow and fluidization: Continuum
and kinetic theory descriptions. Boston, 1994, 211 p.
12. Nilsen M. Modelling of rockfall runout range. Employing
empirical and dynamical methods. Master Thesis in Geo-sciences Discipline: Environmental geology and geo-hazards. Oslo, 2008, 96 p.
22 июня 2016 г.
