CC BY
21
^ Научно-технические ведомости СПбГПУ. Физико-математические науки № 1(165) 2013 ^ УДК 532.59:532.51 7
А.И. Храбрый, Е.М. Смирнов, Д.К. Зайцев
ВЛИЯНИЕ МОДЕЛИ ТУРБУЛЕНТНОСТИ НА РЕЗУЛЬТАТЫ РАСЧЕТА ОБТЕКАНИЯ ПРЕПЯТСТВИЯ ПОТОКОМ ВОДЫ ПОСЛЕ ОБРУШЕНИЯ
ДАМБЫ
A.I. Khrabry, E.M. Smirnov, D.K. Zaytsev
St. Petersburg State Polytechnical University, 29 Politekhnicheskaya St., St. Petersburg, 195251, Russia
INFLUENCE OF TURBULENCE MODEL ON THE RESULTS OF SIMULATION OF DAMBREAK FLOW ABOUT AN OBSTACLE
Представлены результаты численного моделирования обтекания препятствия в форме параллелепипеда потоком, возникшим в результате обрушения дамбы. Показано, что расчеты без использования модели турбулентности переоценивают интенсивность волнообразования. Модель турбулентности k-e обеспечивает более высокий уровень турбулентной вязкости, чем модель SST, и позволяет получить относительно близкое к эксперименту решение.
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ. СВОБОДНАЯ ПОВЕРХНОСТЬ. ОБРУШЕНИЕ ДАМБЫ. ОБТЕКАНИЕ ПРЕПЯТСТВИЯ. ТУРБУЛЕНТНОСТЬ.
The results of numerical simulation of a dambreak flow in the presence of a rectangular obstacle are presented. It is shown that computations without turbulence modeling overestimate the intensity of the wave formation. The k-e turbulence model yields higher turbulent viscosity, than SST model, and provides reasonable agreement with experimental data.
NUMERICAL SIMULATION. FREE SURFACE. DAMBREAK FLOW. OBSTACLE. TURBULENCE
Течения несжимаемой жидкости со свободными поверхностями широко распространены в природе и технике. В последнее время они все более активно исследуются методами вычислительной гидродинамики. Большое внимание, в частности, уделяется изучению течения, возникающего в результате обрушения дамбы. В литературе представлен ряд работ, посвященных численному моделированию течений данного вида в различных условиях: распространение потока в канале постоянного [1] и переменного [2] сечений, натекание потока на вертикальную стенку [3], взаимодействие потока со стоящим на дне канала препятствием в форме клина [4], параллелепипеда [5 — 7] и др. Результаты проведенных расчетов сопоставляются, в той или иной мере, с экспериментальными данными, имеющимися в литературе для соответствующей конфигурации течения.
Вместе с тем, в большинстве публикаций сопоставление с экспериментом проводится без представления сведений, позволяющих оценить качество получаемого численного решения. Вопросам учета эффектов турбулентности в расчетах течений данного вида, равно как и применимости той или иной модели, также не уделяется должного внимания. В то же время есть основания полагать, что учет эффектов турбулентности может оказывать существенное влияние. На это указывает, в частности, обнаруженное в расчетах [8] значительное (свыше 10 %) изменение скорости движения фронта жидкости в зависимости от начальных параметров турбулентности. И это при том, что в статье [8] рассматривался относительно простой случай распространения потока без препятствий; в более сложных течениях степень влияния модели и параметров турбулентности на получаемое решение может оказаться еще
больше.
Настоящая работа посвящена численному моделированию обтекания препятствия в форме параллелепипеда потоком, возникшим в результате обрушения дамбы. Особое внимание уделяется получению решений с пренебрежимо малым влиянием схемных факторов. Результаты расчетов с использованием моделей турбулентности к-г [9] и 88Т [10], а также без учета эффектов турбулентности сравниваются друг с другом и с экспериментальными данными.
Вычислительная методика
Метод VOF. Среди методов, разработанных для решения задач о течении жидкости со свободной поверхностью, наибольшее распространение получил метод VOF (Уо1ите-ОГ-Р1шё) [11], хорошо зарекомендовавший себя как не требовательный к вычислительным ресурсам и способный работать в условиях течений с сильной деформацией свободной поверхности, включая случаи слияния и разделения объемов, занимаемых жидкой и газообразной фазами. Подобное поведение, в частности, имеет место при взаимодействии потока жидкости с плохо обтекаемым препятствием и при опрокидывании потока в результате взаимодействия с вертикальной стенкой.
В данном методе для определения положения межфазной границы, перемещающейся по расчетной сетке, используется распределение специальной «индикаторной» функции С, представляющей собой объемную долю жидкости в ячейках расчетной сетки: при С = 1 ячейка заполнена жидкостью, при С = 0 ячейка пуста. Положению межфазной границы соответствует изоповерхность С = 0,5. Динамика величины С описывается уравнением конвективного переноса:
йС дС дС „
— = — + иI-= 0
йг дг 1 дх1
(1)
где и — компоненты скорости, х1 — пространственные координаты.
В так называемой одножидкостной формулировке метода ^^ используемой в настоящей работе, уравнения движения и переноса параметров турбулентности для системы жидкость-газ записываются как для единой
среды с переменными материальными свойствами, выражаемыми через объемную долю жидкости С:
р=С Рж +(-С )рг; ц = Сцж +(1 -С) '
(2) (3)
где р — плотность среды, ц — динамическая вязкость среды, индексы «ж» и «г» относятся к жидкости и газу соответственно.
Численная схема. В настоящей работе аппроксимация уравнений проводилась по методу конечных объемов, что позволяет использовать неструктурированные расчетные сетки с ячейками различной геометрии. При таком подходе для определения потоков переносимых величин через грани ячейки необходима интерполяция их значений из центров ячеек. Использование стандартных противопоточных аппроксимаций первого или второго порядка для величины С в уравнении (1) приводит к сильному «размытию» межфазной границы из-за численной диффузии. Для уменьшения данного нежелательного эффекта в литературе предложен ряд специальных, так называемых «сжимающих» численных схем. В представленных ниже расчетах для вычисления величины С на грани ячеек использовалась сжимающая схема М-С1С8АМ [13]; для аппроксимации производной по времени использовалась схема Кранка — Николсона. Как показано в работе [12], такая комбинация обеспечивает наименее жесткие ограничения по качеству расчетной сетки и числу Куранта (до 0,75).
Ограничение на шаг по времени для уравнений гидродинамики мягче, чем для уравнения (1). С учетом этого обстоятельства при продвижении по физическому времени использовался алгоритм, предусматривающий несколько дробных шагов для решения уравнения (1) на каждом шаге решения уравнений гидродинамики. Число дробных шагов выбирается таким, чтобы число Куранта для уравнения (1) не превышало 0,55.
Для системы несжимаемых несмешиваемых жидкостей уравнение неразрывности может быть записано в форме условия соленоидаль-ности поля скорости:
ди
дХ;
= 0 .
(4)
Уравнения переноса импульса и параметров турбулентности можно записать в консервативной форме:
д^ + дри = h ,
dt
дх,-
(5)
где ф — переносимая величина (компонента скорости и, или параметр турбулентности), ЯИ5 — правая часть уравнения, которая может включать диффузионные слагаемые и источни-ковые члены.
Данная форма может использоваться только в условиях согласованности аппроксимаций для плотности среды и объемной доли Сна гранях ячеек. Поскольку на границе жидкость-газ плотность терпит разрыв, изменяясь на несколько порядков, малейшая несогласованность аппроксимаций приводит к возникновению больших погрешностей в уравнении (5), что, в частности, может приводить к искажению формы свободной поверхности. В связи с этим была предложена иная форма записи, оказавшаяся менее чувствительной к способу аппроксимации плотности на грань:
рЗф + дрфщ _фдрщ = н (6)
дt дх^ дхг-
Для аппроксимации по времени использовалась схема Кранка — Николсона (7), приводящая к существенно меньшим искажениям свободной поверхности, чем неявная трехслойная схема:
рфп _ф"_1 + &_ф= яжГф). (7)
Дt дх^ дх.1
Здесь индексы п — 1 и п относятся к старому и новому слоям во времени соответственно. Величины р , ф , и относятся к промежуточному слою ^_1/2. Значения плотности на этом слое вычисляются как полусумма значений на слоях п — 1 и п:
рп _1 +рп Р 2 '
В то же время для вычисления ф и и подобный метод, как оказалось, не подходит, так как с ним схема становится неустойчивой во времени. Для преодоления неустойчивости в выражения для ф и и вводились стабилизирующие поправки второго порядка малости:
Ф = Ф" +2Ф"-1 + Аф, Аф = ( - 2 Фи-! + Фи )16
и аналогично для й.
Постановка задачи
Проводится численное моделирование обтекания препятствия в форме параллелепипеда потоком, возникшим в результате обрушения дамбы. Постановка задачи соответствует эксперименту, проведенному в Морском научно-исследовательском институте Нидерландов (MARIN); результаты измерений приводятся, в частности, в работе [7]. В эксперименте использовался бак в форме параллелепипеда с основанием 3,2 х 1,0 м (рис. 1). На дне резервуара, на расстоянии 0,68 м от правой на рисунке торцевой стенки, располагалось препятствие размерами 0,16 х 0,16 х 0,40 м, на равном расстоянии от боковых стенок. В начальный момент времени вода занимает объем 1,2 х 0,55 х 1,0 м у левой торцевой стенки (в эксперименте воду удерживала перегородка, которая в начальный момент времени резко убиралась). Поток распространяется по баку, взаимодействуя с препятствием, отражается от торцевой стенки (в качестве иллюстрации на рис. 1 приведено расчетное положение свободной поверхности на момент времени 0,64 с), затем вторично проходит по резервуару и отражается от левой торцевой стенки. Проводились измерения высоты уровня жидкости во времени в точке H2, находящейся в плоскости симметрии, на расстоянии 1 м от правой торцевой стенки.
Рис. 1. Схема эксперимента MARIN и расчетное положение свободной поверхности на момент времени 0,64 с;
H2 — точка мониторинга высоты жидкости. Все размеры даны в миллиметрах
Результаты расчетов
Выбор расчетной сетки. Для выбора расчетной сетки и шага по времени предварительно решалась двумерная задача, соответствующая плоскости симметрии исходной задачи (без учета эффектов турбулентности). Расчеты проводились на трех декартовых сетках размерностями 160 х 50, 320 х 100 и 640 х 200 ячеек (далее по тексту употреблены названия грубая, средняя и мелкая сетки). Шаг по времени в каждом из расчетов был постоянным (0,004, 0,002 и 0,001 с соответственно). При таких шагах типичное число Куранта в каждом из расчетов составляло 0,8 (в отдельные моменты времени локальные значения числа Куранта достигали нескольких единиц). Также для средней сетки был проведен расчет с шагом по времени, уменьшенным в два раза (0,001 с).
Расчетная форма объема, занимаемого жидкостью вблизи препятствия по прошествии 1,2 с, для данных четырех расчетов приведена на рис. 2. Можно видеть, что решения, полученные на средней и мелкой сетках, достаточно близки, а решение на грубой сетке заметно от них отличается. Таким образом, решение, полученное на средней сетке, можно считать сеточ-но-сошедшимся. Также видно, что результаты очень близки при использовании стандартного и уменьшенного шага по времени. Таким образом, шаг по времени 0,002 с достаточен для проведения расчетов на средней сетке.
Влияние модели турбулентности. Двумерные расчеты. Для оценки характера и степени вли-
2,2 2,4 2,5 2.В 3.0 я,м
Рис. 2. Положение свободной поверхности на момент времени 1,2 с (2D-расчет): средняя (1, 2) расчетная сетка со стандартным (1) и уменьшенным (2) шагами по времени; грубая (3) и мелкая (4) сетки
яния модели турбулентности на исследуемое течение были проведены двумерные расчеты с использованием моделей турбулентности k- s [9] и SST [10], а также без учета турбулентности. Результаты для момента времени 1,44 с приведены на рис. 3,а. Видно, что при расчете, не учитывающем эффектов турбулентности, имеет место интенсивное образование волновых структур на поверхности жидкости. При этом из-за опрокидывания волн высота жидкости в точке H2 может определяться неоднозначно: высотой можно считать как верхнее, так и нижнее положение свободной поверхности над точкой H2. На рис. 3,б приведены зависимости от времени для «верхней» и «нижней» высот в расчете без учета турбулентности.
Поскольку в действительности гребни опрокидывающихся волн, возникающих при ударе о препятствие или стенку, разлетаются на брызги, нижнее положение свободной поверхности представляется более правильным. В связи с этим на последующих графиках приведены высоты, замеренные по нижнему положению свободной поверхности.
На рис. 4 приведены зависимости высоты жидкости в точке H2 от времени в разных вариантах расчетов. Видно, что без учета турбулентности и с моделью SST течение за весь период времени получается значительно более хаотичным, чем при расчете с моделью k- s . Таким образом, модель турбулентности k- s обеспечивает достаточно высокий уровень турбулентной вязкости для подавления и сглаживания волн, «гуляющих» по расчетной области. Модель SST дает значительно более низкий уровень турбулентной вязкости, чем модель ks, и решение получается близким к расчету без учета турбулентности.
Влияние модели турбулентности. Трехмерные расчеты. С учетом симметрии задачи расчетная область представляла собой половину экспериментальной. Расчетная сетка имела 320 х 100 ячеек вдоль плоскости симметрии и 35 ячеек в поперечном направлении.
С учетом выводов из решения двумерных задач расчеты в полной трехмерной постановке, соответствующей эксперименту MARIN, проводились в двух «крайних» постановках: без учета эффектов турбулентности и с использованием модели k- s . Результаты представлены
¡ш
2
2.2
Рис. 3. Результаты двумерных расчетов: а — положение жидкости в момент времени 1,44 с без учета (1) и с учетом (2) турбулентности (модель к- е); б — зависимость от времени высоты жидкости в точке Н2 по верхнему (3) и нижнему (4) положениям свободной поверхности (без учета турбулентности)
Рис. 4. Зависимости от времени высоты жидкости в точке Н2, полученные без учета турбулентности (1) и с использованием моделей SST (2) и к- е (3)
на рис. 5. Видно, что при использовании модели к- е достигается существенно лучшее согласование с экспериментальными данными [7]. Исключение из вычислительной модели эффектов турбулентности приводит к чрезмерно интенсивному волнообразованию. Примечательно также, что расчет с моделью к- е правильнее предсказывает подъем уровня жидкости на начальной стадии взаимодействия потока с препятствием (? < 1,5 с) и время прихода волны, отразившейся от левой стенки бака (?« 4,7 с).
Рис. 5. Зависимости от времени высоты жидкости в точке H2, полученные по эксперименту MARIN (1), 3D-расчету без учета турбулентности (2) и с моделью к- е (3)
По результатам проведенного исследования можно заключить, что учет турбулентности может оказывать существенное влияние на получаемое решение в сложных течениях жидкости со свободной поверхностью. При этом модель турбулентности k- е обеспечивает более высокий уровень турбулентной вязкости, чем модель SST, и позволяет получить решение, относительно близкое к эксперименту. Расчеты без использования модели турбулентности переоценивают интенсивность волнообразования.
СПИСОК ЛИТЕРАТУРЫ
1. Hogg, P.W. An implicit algorithm for capturing method [Текст] / P.W. Hogg, X.J. Gu, D.R.Emerson // sharp fluid interfaces in the volume of fluid advection Technical Report DL-TR-2006-001. — 2006. — 26 p.
2. Kocaman, S. The effect of lateral channel contraction on dam break flows: Laboratory experiment [Текст] / S. Kocaman, H. Ozmen-Cagatay // Journal of Hydrology. - 2012. - Vol. 432-433. - P. 145-153.
3. Park, I.R. A volume-of-fluid method for incompressible free surface flows [Текст] / I.R. Park, K.S. Kim, J. Kim, S.H. Van // Int. J. Numer. Meth. Fluids. - 2009. -Vol. 61. - P. 1331-1362.
4. Biscarini, C. CFD modelling approach for dam breakflow studies [Текст] / C. Biscarini, S. Di Francesco, P. Manciola // Hydrology and Earth System Sciences. -2010. - Vol. 14. - P. 705-718.
5. Минаков, А.В. Численное моделирование течений вязкой несжимаемой жидкости с подвижными границами [Текст]: Дис. ... канд. физ.-мат. наук: 05.13.18. Защищена 19.12.08/ Минаков Андрей Викторович. - Красноярск, Сибирский федеральный университет, 2008. - 189 с. - Библиогр.: с. 175-189.
6. Wemmenhove, R. Numerical simulation of two-phase flow in offshore environments [Текст]: PhD thesis/ Wemmenhove Rik. — University of Groningen, 2008. -141 p. - Bibliogr.: P. 121-125.
7. Kleefsman, K.M.T. A Volume-of-Fluid based simulation method for wave impact problems [Текст] / K.M.T. Kleefsman, G. Fekken, A.E.P. Veldman // Journal of Computational Physics. - 2005. - Vol. 206. - P. 363393.
8. Park, I.R. Numerical investigation of the effects of turbulence intensity on dam-break flows [Текст] / I.R. Park, K.S. Kim, J. Kim, S.H. Van // Ocean Engineering. - 2012. - Vol. 42. - P. 176-187.
9. Launder, B.E. The numerical computation of turbulent flows [Текст] / B.E. Launder, D.B. Spalding // Computer Methods in Applied Mechanics and Engineering. - 1974. - Vol. 3. - P. 269-289.
10. Menter, F.R. Ten years of industrial experience with the SST turbulence model [Электронный ресурс] / F.R. Menter, R. Langtry, M. Kuntz // Turbulence, Heat and Mass Transfer IV, Antalya, Turkey, October 12— 17, 2003. - CD-ROM proceedings. - Begell House, Inc., 2003. - 8 p.
11. Hirt, C.W. Volume of fluid (VOF). Method for the dynamics of free boundaries [Текст] / C.W. Hirt, B.D. Nichols // Journal of Computational Physics. -1981. - Vol. 39. - P. 201-226.
12. Khrabry, A.I. Solving the convective transport equation with several high-resolution finite volume schemes: Test computations [Текст] / A.I. Khrabry, E.M. Smirnov, D.K. Zaytsev // Computational Fluid Dynamics 2010. -New-York: Springer, 2011. - P. 535—540.
13. Wactawczyk, T. Remarks on prediction of wave drag using VOF method with interface capturing approach [Текст] / T. Wacfowczyk, T Koronowicz // Archives ofCivil and Mechanical Engineering. - 2008. - Vol. 8. - P. 5—14.
ХРАБРЫЙ Александр Иосифович — ведущий программист кафедры гидроаэродинамики Санкт-Петербургского государственного политехнического университета.
195251, г. Санкт-Петербург, Политехническая ул., 29 xbr@list.ru
СМИРНОВ Евгений Михайлович — доктор физико-математических наук, профессор, заведующий кафедрой гидроаэродинамики Санкт-Петербургского государственного политехнического университета, 195251, г. Санкт-Петербург, Политехническая ул., 29 aero@phmf.spbstu.ru
ЗАЙЦЕВ Дмитрий Кириллович — кандидат физико-математических наук, доцент кафедры гидроаэродинамики Санкт-Петербургского государственного политехнического университета. 195251, г. Санкт-Петербург, Политехническая ул., 29
© Санкт-Петербургский государственный политехнический университет, 2013
