ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЛАВОВЫХ ПОТОКОВ В МОДЕЛЯХ ИЗОТЕРМАЛЬНОЙ ВЯЗКОЙ МНОГОФАЗНОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЛАВОВЫХ ПОТОКОВ В МОДЕЛЯХ ИЗОТЕРМАЛЬНОЙ ВЯЗКОЙ МНОГОФАЗНОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ
Научная статья
1, 2 Институт математики и механики им. Н.Н. Красовского уральского отделения РАН, Екатеринбург, Россия
* Корреспондирующий автор (starodubtsevayv[at]ya.ru)
АннотацияПоток лавы начинает формироваться, когда расплавленная порода извергается на поверхность Земли и медленно распространяется на поверхности из вулканического жерла. Извержения создают различные потоки лавы (например, течения различной структуры и скорости потока) под действием гравитации в зависимости от химического состава, температуры магматических пород и топографии поверхности, по которой протекает лава. Несмотря на то, что потоки вулканической лавы не оказывают существенного влияния на жизнь людей, их опасность немалая, поскольку горячая лава убивает растительность, разрушает инфраструктуру и может вызвать наводнение из-за таяния снега/льда. Следуя развитию вычислительных ресурсов, численное моделирование лавовых потоков за последние несколько десятилетий продвигалось от моделирования одномерных потоков к моделированию трехмерных потоков, которое наиболее адекватно способно отразить реальные природные процессы. В этой статье разрабатываются трехмерные численные модели течений изотермальной вязкой ньютоновской многофазной жидкости на различных поверхностях под действием силы тяжести, с целью исследовать динамику и взаимодействие потоков лавы. Полное моделирование потока лавы является сложной задачей с физической, математической и численной точек зрения. Математическая модель включает в себя уравнение Навье—Стокса, уравнение несжимаемости и уравнения переноса фаз с соответствующими начальными и граничными условиями. Для численной аппроксимации математической модели применяется метод конечных объемов. Программные коды реализованы в пакете ANSYS Fluent на языке C. При проведении численных экспериментов использовалась ЭВМ параллельного действия. Демонстрируются результаты расчетов модельного эксперимента. Модели реконструкции потоков лавы могут оказывать существенную помощь при проектировании барьеров, отражающих потоки лав. Доступность технологических и научных данных (например, данные спутниковых мониторингов, высокоскоростные алгоритмы расчетов и реалистичные модели) позволят интегрировать данные в модели с традиционными методами изучения вулканической активности, что позволит более эффективно использовать результаты.
Ключевые слова: вязкая жидкость, многофазная жидкость, уравнение Навье–Стокса, краевая задача, численное моделирование.
NUMERICAL SIMULATION OF LAVA FLOWS IN MODELS OF ISOTHERMAL VISCOUS MULTIPHASE INCOMPRESSIBLE FLUID
Research article
Korotkiy A.I.1, Tsepelev I.A.2
1, 2 N.N. Krasovsky Institute of Mathematics and Mechanics of the Ural Branch of the Russian Academy of Sciences, Yekaterinburg, Russia
* Corresponding author (starodubtsevayv[at]ya.ru)
AbstractA lava flow begins to form when molten rock erupts onto the surface of the Earth and slowly spreads to the surface from a fissure vent. Eruptions create different lava flows (for example, flows of a different structure and flow velocity) under the influence of gravity, depending on the chemical composition, temperature of igneous rocks, and the topography of the surface over which the lava flows. Despite the fact that volcanic lava flows do not have a significant impact on people's lives, their danger is considerable, since hot lava kills vegetation, destroys infrastructure, and can cause flooding due to melting snow/ice. Following the development of computing resources, numerical modeling of lava flows over the past few decades has moved from modeling one-dimensional flows to modeling three-dimensional flows, which is most adequately able to reflect real natural processes. In order to investigate the dynamics and interaction of lava flows, the current article develops three-dimensional numerical models of flows of an isothermal viscous Newtonian multiphase fluid on various surfaces under the influence of gravity. A complete simulation of a lava flow is a challenging task from a physical, mathematical, and numerical point of view. The mathematical model includes the Navier-Stokes equation, the incompressibility equation, and the phase transfer equations with corresponding initial and boundary conditions. The finite volume method is used for numerical approximation of the mathematical model. The program codes are implemented in the ANSYS Fluent package in C. When conducting numerical experiments, a parallel-acting computer was used. The article demonstrates the results of calculations of the model experiment. Lava flow reconstruction models can provide significant assistance in the design of barriers reflecting lava flows. The availability of technological and scientific data (such as satellite monitoring data, high-speed calculation algorithms, and realistic models) will allow for integrating data into models with traditional methods of studying volcanic activity, which will allow more efficient use of the results.
Keywords: viscous liquid, multiphase liquid, Navier-Stokes equation, boundary value problem, numerical modeling.
ВведениеЛавовые потоки – это потоки частично расплавленной породы, которые распространяются по поверхности Земли под действием силы тяжести. При охлаждении такого расплава на его поверхности образуется твердая корка, расплав постепенно затвердевает, пока не превратится в твердое тело. Математическое и компьютерное моделирование играют важную роль в понимании структур потока лавы, его морфологии и тепловой эволюции (например, см., [1], [2]). Следуя развитию вычислительных ресурсов, численное моделирование лавовых потоков за последние несколько десятилетий продвигалось от моделирования одномерных потоков к моделированию трехмерных потоков, которое наиболее адекватно способно отразить реальные природные процессы. В работе [2] представлен обзор численных подходов и программных средств для моделирования распределения потоков лавы. С развитием средств космического и аэромониторинга, в моделях реконструкции потоков лавы могут использоваться цифровые данных об изменении поверхности местности под влиянием накопления лавы, и тем самым оказывать помощь при проектировании барьеров, отражающих лавовую опасность. Доступность данных спутникового мониторинга, высокоскоростных ЭВМ, реалистичных математических моделей, позволяют интегрировать данные и модели, и эффективно использовать полученные результаты. В данной работе для моделирования движения лавы используется метод переноса объёмной доли (VOF). Этот метод основан на дискретизации расчетной области неструктурированной сеткой. Для аппроксимации систем дифференциальных уравнений, которые описывают движение ньютоновской вязкой несжимаемой жидкости под действием внешних сил [3], [4] и учитывает влияние теплопереноса на границе взаимодействия различных жидкостей, применяется метод конечных объемов [5].
Реология лавы зависит от ее химического состава, температуры, содержания кристаллов и водяного пара. Так, например, лавы с преобладанием базальта более плотные и вязкие, чем лавы с преобладанием кремния. Реология зависит также от времени как результат охлаждения, кристаллизации и конденсации водяного пара. Охлаждение определяется сочетанием переноса тепла в жидкой лаве с переносом его на поверхность лавы и в окружающую атмосферу посредством радиационного и конвективного теплообмена [6]. В работе [7] численное моделирование течения лавовых потоков осуществлялось в двумерной постановке в рамках движения теплопроводной двухфазной неньютоновской несжимаемой жидкости с учетом нелинейного теплообмена на границе жидкостьвоздух.
В этом исследовании представлена трехмерная модель потока вязкой несжимаемой жидкости, находящейся под действие силы тяжести, которая может описывать изотермические течения лавы. В математической модели течения лавовых потоков, представленной в настоящем исследовании, не учитываются явно тепловые и реологические процессы, которые, несомненно, влияют на динамику распространения лавы. Также модель не учитывает процессы дегазации и зависимость плотности от температуры. Вместо этого используется некоторая модельная зависимость физических характеристик лавового потока (а именно, плотность и вязкость) от времени, чтобы симулировать процесс затвердевания расплавленных пород.
Постановка задачи
Многофазный поток – это поток с участием более чем одной жидкости. Каждая такая жидкость обладает своими физическими характеристиками, а именно, плотностью и вязкостью. Жидкости не смешиваются, но обладают некоторой областью межфазного взаимодействия и могут взаимодействовать друг с другом посредством межфазных сил, например, силы поверхностного натяжения в зоне взаимодействия.
Опишем математическую модель рассматриваемого движения жидкости. В качестве основных уравнений состояния жидкости примем уравнения движения двухфазной вязкой жидкости. В модельной области (см. рисунок 1) движение такой вязкой жидкости на промежутке времени – начало времени наблюдения за процессом, – конечный момент времени наблюдения за процессом, представляется уравнением Навье–Стокса, и уравнением несжимаемости (неразрывности) [3], [4]
плотность и вязкость среды может быть вычислена из соотношенийЗдесь – вязкость и плотность жидкостей, которые поступают в рассматриваемую область извне; – вязкость и плотность жидкости, которая заполняют модельную область в начальный момент времени. Модельная область ограничена снизу поверхностью исследуемой местности – горизонтальные размеры региона и – высота модельной области. На нижней границе модельной границе и боковых границах задается условие прилипания: . На той части нижней поверхности, где расположено жерло, задаются условия: где – скорость экструзии магмы. На верхней границе модельной области задаются условия . Заметим, что указанные граничные условия гарантируют выполнение условия (2).
Задача состоит в нахождении функций
являющихся решением начально-краевой задачи (1) -- (4) с указанными выше начальными и граничными условиями. Величины считаются известными. Данная задача нелинейна относительно искомых функций, аналитическое нахождение решений такой задачи математически сложно, поэтому альтернативным методом исследования является численное моделирование. В работе [8] проводится исследование корректности постановок начально-краевых задач для упомянутых выше моделей. Вводится понятие обобщенного (слабого) решения соответствующей математической задачи. Также исследуется разрешимость задачи и указана зависимость решения задачи от исходных данных и параметров модели.
Система уравнений (1)–(6) имеет устойчивое решение, если произвольное по отношению к нему малое возмущение, возникнув в какой-либо момент, далее со временем убывает. В противном случае решение неустойчиво. Можно показать, что при малых числах Рейнольдса, когда вязкие силы играют преобладающую роль, гидродинамические решения устойчивы. Неустойчивости возникают при больших числах Рейнольдса, когда существенны инерционные силы. Многофазная жидкость также будет находиться в неустойчивом состоянии, если сверху находятся более плотные фазы жидкости. Такая система будет эволюционировать под действием силы тяготения. В данном исследовании наверху находится фаза, моделирующая воздух, что не позволяет силе тяжести вывести систему из устойчивого равновесия.
Рис. 1 – Модельная область и ее граница
Численный метод решения задачиДля численного анализа модельной задачи мы используем программное обеспечение Ansys Fluent (https://www.ansys.com/products/fluids/ansys-fluent), основанное на методе конечных объемов. Модельная область разбивается на гексаэдров, составляющих конечные объемы. Численные коды используют многофазную нестационарную VOF (Volume Of Fluid) модель [9], а решатель использует неявную по времени схему интегрирования уравнений (1) -- (4) для совместного определения поля скоростей, давления и объемной доли жидкости. Метод VOF основан на технике отслеживания поверхности, применяемой к фиксированной эйлеровой сетке. Он разработан для двух или более несмешивающихся жидкостей, когда положение границы раздела между жидкостями четко определяется. Уравнения импульса (1) являются общими для всех жидкостей, и объемная доля каждой из жидкостей (4) в каждой расчетной ячейке (элемент конечного объема) отслеживается по всей области. Если в i ячейке расчетной сетки , тогда эта ячейка полностью заполнена p жидкостью. Если , тогда в ячейке i жидкость p не присутствует. Если , то ячейка i есть зона интерфейса между жидкостью p и другими жидкостями. Для всех моментов времени для каждой ячейки сетки выполнено условие (6).
Для аппроксимации давления и лапласиана используются численные схемы второго порядка точности; для дискретизации конвективных членов -- монотонные схемы (см., например, [10]). Давление дискретизируется по ступенчатой схеме второго порядка PRESTO! (PREssure STaggering Option) [11]. Давление и скорость вычислялись SIMPLE (Semi-Implicit Method for Pressure Linked Equations) методом [12], где параметры релаксации выбраны равными 0.01 и 0.3 для скорости и давления, соответственно. Учитывая нелинейность задачи, временной шаг выбирается в диапазоне от 1с. до 10с. в зависимости от устойчивости и для оптимизации скорости сходимости систем линейных алгебраических уравнений (СЛАУ), полученных после дискретизации задач (1), (2), (4). Неявная схема позволяет выполнять устойчивые вычисления с относительно большим временным шагом. Для решения всех СЛАУ используется многосеточные методы [10].
Заметим, что число Рейнольдса ([13], C.87) для лавы мало, поскольку лава движется медленно и вязкость ее большая. Поэтому течение лавы всегда ламинарное. При этом число Рейнольдса для воздуха большое. Отношение вязкости лавы к вязкости воздуха может достигать порядка . Численная неустойчивость, которая может возникать в процессе вычислений, приводит к тому, что движение в воздухе становиться турбулентным, что делает дальнейшее проведение расчетов невозможным. Для преодоления этой трудности приходится уменьшать параметр релаксации для определения скорости и проводить расчеты с очень малым шагом по времени. Отметим, что попытки решать данную задачу явными методами не увенчались успехом.
Учитывая высокую размерность дискретной задачи, для вычислений привлекался вычислительный кластер “Уран” (ИММ УрО РАН, г.Екатеринбург). Вычисления проводились на 1 узле, который состоит их двух 18-и ядерных процессоров Intel(R) Xeon(R) CPU E5-2697 v4 2.30GHz; оперативная память 256 GB; кэш-память 45 MB SmartCache. Ориентировочное время счета составило около 36 часов на 32 ядрах CPU.
Вычислительный экспериментРасчетная сетка строится по следующему алгоритму. В прямоугольнике строится равномерная прямоугольная сетка . Далее определяется гексаэдры (блоки), каждый из которых определяется 8 вершинами:
Каждый такой блок разбивается в вертикальном направлении точками таким образом, что высота каждой ячейки определяется по формуле
Такое разбиение делает высоту ячеек, которые ближе к нижней границе, кратно меньше, чем у тех ячеек, которые ближе к верхней границе. Таким образом, удается аппроксимировать лавовые потоки более качественно учетом того, что высота лавы уменьшается по мере изменения рельефа местности весьма существенно. Для генерации сетки использовался пакет OpenFOAM (https://www.openfoam.com/). Число неортогональности сетки в данном расчет не превосходило 17, а среднее число неортогональности равнялось 3. Построенная таким образом гексаэдральная сетка может быть легко масштабирована в любом направлении практически без потери качества. При уменьшении вязкости жидкости, высота растекания жидкости уменьшается. В этом случае просто уменьшаются вертикальные высоты ячеек сетки, примыкающие к нижней поверхности, для более качественного моделирования тонких структур.
Рис. 2 – Ортогональность – угловое отклонение вектора (в градусах), соединяющего геометрические центры двух ячеек, и вектора нормали к их общей грани. Иногда под ортогональностью понимают значение косинуса этого угла
Рис. 3 – Построение расчетной сетки
Примечание: здесь
До начала эксперимента предполагается, что модельная область заполнена жидкостью с характеристиками оторая моделирует воздушную среду. С момента времени жидкость с характеристиками (см. рисунок 1) из некоторой области нижней границы втекает в область модели со скоростью -- вектор внешней нормали в соответствующей точке нижней границы. Жидкость стекает по рельефу в течение часов, после чего подача жидкости 1 прекращается (см. рисунок 4). С момента времени мы изменяем свойства жидкости 1, увеличивая вязкость излившейся жидкости в 10 раз, После этого мы продолжаем вычисления, и с момента времени жидкость с характеристиками из той же области на нижней границе и с той же скоростью начинает поступать внутрь модельной области. Жидкость 2 течет над жидкостью 1 в течение примерно часов, после чего подача жидкости 2 также прекращается (см. рисунок 5). Жидкость 2 практически полностью перекрывает жидкость 1, поскольку вязкость жидкости 1 выше и она движется медленнее жидкости 2 в результате действия силы тяжести и изостатического сжатия жидкости 1 жидкостью 2. С момента времени мы изменяем свойства жидкости 1 и 2, увеличивая вязкость излившейся жидкости в 10 раз, Вычисления продолжаются и с момента времени жидкость с характеристиками из той же области на нижней границе и со скоростью поступает внутрь модельной области. В начале жидкость 3 находится под жидкостями 1 и 2. Далее жидкость 3 выходит на поверхность и течет поверх двух более «старых» жидкостей, толкая их вниз и еще больше сжимая нижележащие жидкости. Вычисления продолжаются до момента времени часов (см. рисунок 6). Такой сценарий излияния расплавленных пород генерирует слоистый по своему составу поток лавы.
Рис. 4 – Демонстрируются границы раздела между воздухом и жидкостью 1, характеризующие движение этой жидкости по рельефу в моменты времени 0s, 3000s, 6000s и 9000s
Рис. 5 – Совместное движение жидкости 1 и жидкости 2 по рельефу в моменты времени 10000s, 18000s, 26000s и 33000s
Рис. 6 – Совместное движение жидкости 1, жидкости 2 и жидкости 3 по склону в моменты времени 34000s, 43000s, 53000s и 63000s
Таблица 1 – Физические параметры жидкостей в численном эксперименте
ЗаключениеВ данной работе представлена численная реализация пространственной мультифазной модели движения ньютоновских вязких несжимаемых жидкостей, которые распространяются в воздушной среде под действие силы тяжести. Проведены расчеты модельного примера движения таких жидкостей, с учетом того, что каждая жидкость имеет свои физические характеристики. Учитывая применение пакета ANSYS Fluent для реализации решателей систем дифференциальных уравнений, данный подход может быть применен для моделирования более сложных задач в пространственных областях. Моделирование, представленное в работе, способно предсказать динамику и морфологию лавового потока, что в сочетании уже с традиционными методами оценки лавовой опасности, поможет уменьшить риски экономических потерь при извержении вулканов, особенно в густонаселенных регионах мира ([14]).
Финансирование Работа выполнена при поддержке РФФИ и DFG в соответствии с исследовательским проектом №205112002. При проведении работ был использован суперкомпьютер «Уран» ИММ УрО РАН. | Funding This work was supported by RFBR and DFG according to the research project No 205112002. Our work was performed using ”Uran” supercomputer of IMM UB RAS. |
Конфликт интересов Не указан. | Conflict of Interest None declared. |
Список литературы / References
- Macedonio G. A computer model for volcanic ash fallout and assessment of subsequent hazard / G. Macedonio, A. Costa, A. Longo // Computers & Geosciences, 2005. — Vol. 31. — P. 837–845. DOI: 10.1016/j.cageo.2005.01.013 .
- Cordonnier B. Benchmarking lavaflow models / B. Cordonnier, E. Lev, F. Garel // Geological Society, London, Special Publications, 2015. — Vol. 426. — P. 425–445. DOI: 10.1144/SP426.7A .
- Kolev N.I. Multiphase flow dynamics/ N.I. Kolev // Berlin, Heidelberg: SpringerVerlag, 2011. — 781 p.
- Нигматулин Р.И. Динамика многофазных сред. Часть 1 / Р.И. Нигматулин. — М.: Наука, 1987. — 464 с.
- Ferziger J.H. Computational Methods for Fluid Dynamics / J.M. Ferziger, M. Perit // SpringerVerlag Berlin Heidelberg, 2002. — 426 p. DOI: 10.1007/9783642560262 .
- Griffiths W.R. The Dynamics of Lava Flows / W.R. Griffiths // Annual Review of Fluid Mechanics, 2000. — Vol. 32.— P. 477–518. DOI: 10.1146/annurev.fluid.32.1.477 .
- Tsepelev I. Computational modelling of lava flows in multiphase fluid models with temperature and shear ratedependent rheology / I. Tsepelev, A. Korotkii // AIP Conference Proceedings, 2020. —Vol. 2312, No. 050022. — 8 p. DOI: https://doi.org/10.1063/5.0035525 .
- Цепелев И.А. Гравитационное течение двухфазной вязкой несжимаемой жидкости / И.А. Цепелев, А.И. Короткий, Ю.В. Стародубцева // Труды ИММ УрО РАН, — 2021, — Том 27, N
- Hirt C.W. Volume of fluid (VOF) method for the dynamics of free boundaries / C.W. Hirt, B.D. Nichols // J. Comput.Phys, 1981. — Vol. 39, No. 1, — P. 201–225.
- IsmailZadeh A. Computational Methods for Geodynamics / A. IsmailZadeh, P. Tackley. — Cambridge University Press, Cambridge, 2010. — 347 p. DOI: 10.1017/CBO9780511780820 .
- Peyret R. Handbook of Computational Fluid Mechanics / R. Peyret. — Academic Press Limited, USA, 1996. —467 p.
- Patankar S.V. A calculation procedure for heat and mass transfer in threedimensional parabolic flows / S.V. Patankar, D.B. Spalding // Int. J. Heat Mass Transfer, 1972. — Vol. 15, — P. 1787–1806. DOI: 10.1016/00179310(72)900543 .
- Chandrasekhar S. Hydrodynamic and hydromagnetic stability / S. Chandrasekhar. –– Oxford: Clarendon Press, 1961. ––654 p.
- IsmailZadeh A. Integrating natural hazard science with disaster risk reduction policy / A. IsmailZadeh // Advancing Culture of Living with Landslides”, 2017. — Springer, Cham, — P. 167–172. DOI: 0.1007/978331959469913.
Список литературы на английском языке / References in English
- Macedonio G. A computer model for volcanic ash fallout and assessment of subsequent hazard / G. Macedonio, A. Costa, A. Longo // Computers & Geosciences, 2005. — Vol. 31. — P. 837–845. DOI: 10.1016/j.cageo.2005.01.013 .
- Cordonnier B. Benchmarking lavaflow models / B. Cordonnier, E. Lev, F. Garel // Geological Society, London, Special Publications, 2015. — Vol. 426. — P. 425–445. DOI: 10.1144/SP426.7A .
- Kolev N.I. Multiphase flow dynamics/ N.I. Kolev // Berlin, Heidelberg: SpringerVerlag, 2011. — 781 p.
- Nigmatulin R.I. Dinamika mnogofaznyh sred [Dynamics of multiphase media]. Part 1 / R.I. Nigmatulin. - M.: Nauka, 1987— - 464 p. [in Russian]
- Ferziger J.H. Computational Methods for Fluid Dynamics / J.M. Ferziger, M. Perit // SpringerVerlag Berlin Heidelberg, 2002. — 426 p. DOI: 10.1007/9783642560262 .
- Griffiths W.R. The Dynamics of Lava Flows / W.R. Griffiths // Annual Review of Fluid Mechanics, 2000. — Vol. 32.— P. 477–518. DOI: 10.1146/annurev.fluid.32.1.477 .
- Tsepelev I. Computational modelling of lava flows in multiphase fluid models with temperature and shear ratedependent rheology / I. Tsepelev, A. Korotkii // AIP Conference Proceedings, 2020. —Vol. 2312, No. 050022. — 8 p. DOI: https://doi.org/10.1063/5.0035525 .
- Tsepelev I.A. Gravitacionnoe techenie dvuhfaznoj vjazkoj neszhimaemoj zhidkosti [Gravitational flow of a two-phase viscous incompressible fluid] / I.A. Tsepelev, A.I. Korotkiy, Yu.V. Starodubtseva // Trudy IMM UrO RAN [Proceedings of IMM UrO RAS], - 2021, - Volume 27, N 4. [in Russian]
- Hirt C.W. Volume of fluid (VOF) method for the dynamics of free boundaries / C.W. Hirt, B.D. Nichols // J. Comput.Phys, 1981. — Vol. 39, No. 1, — P. 201–225.
- IsmailZadeh A. Computational Methods for Geodynamics / A. IsmailZadeh, P. Tackley. — Cambridge University Press, Cambridge, 2010. — 347 p. DOI: 10.1017/CBO9780511780820.
- Peyret R. Handbook of Computational Fluid Mechanics / R. Peyret. — Academic Press Limited, USA, 1996. —467 p.
- Patankar S.V. A calculation procedure for heat and mass transfer in threedimensional parabolic flows / S.V. Patankar, D.B. Spalding // Int. J. Heat Mass Transfer, 1972. — Vol. 15, — P. 1787–1806. DOI: 10.1016/00179310(72)900543.
- Chandrasekhar S. Hydrodynamic and hydromagnetic stability / S. Chandrasekhar. –– Oxford: Clarendon Press, 1961. ––654 p.
- IsmailZadeh A. Integrating natural hazard science with disaster risk reduction policy / A. IsmailZadeh // Advancing Culture of Living with Landslides”, 2017. — Springer, Cham, — P. 167–172. DOI: 0.1007/978331959469913.