CC BY
29
УДК 502/5045:627.13:519.852.6
М. А. ВОЛЫНОВ, И. В. ГУГУШВИЛИ
Государственное научное учреждение
Всероссийский научно-исследовательский институт гидротехники и мелиорации имени А. Н. Костикова
ОБ АДЕКВАТНОСТИ МЕТОДА 3D ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ВОЛНЫ ПРОРЫВА
В статье рассматривается задача 3D численного моделирования обрушения столба жидкости и распространения волны прорыва по мокрому руслу. Для моделирования использованы полные уравнения Навье-Стокса. Турбулентность описывается моделью динамики больших вихрей с консервативным методом трассирования свободной поверхности. Возможность решения такой задачи за приемлемое время обусловлена применением технологии использования графических вычислителей. Результаты численного моделирования сопоставляются с экспериментальными данными.
Волны прорыва, 3D моделирование, уравнения Навье-Стокса, физический эксперимент, технология CUDA, модель динамики больших вихрей, расчетная сетка, свободная поверхность.
The article considers a task of the 3D numerical simulation of the liquid column failure and break wave propagation on the sea wetted bed. Complete Navier-Stokes equations were used for simulation. Turbulence is described by a model of big vortexes dynamics with a conservative method of free surface tracing. The possibility of solving such a task during the acceptable time is conditioned by a usage the graphic calculators technologies. The results of numerical simulation are compared with experimental data.
Break waves, 3D simulation, Navier-Stokes equations, physical experiment, CUDA technology, model of big vortexes dynamics, computational grid, free surface.
Прогнозирование последствий проры- параметров волн прорыва в непосредствен-
ва напорного фронта гидроузла, в частности расчет параметров волны прорыва и зоны возможного затопления, является рутинным требованием при проектировании и эксплуатации ГТС. На сегодняшний день для такого прогноза существует большое количество программных комплексов (ПК), основанных как на одномерной, так и на двумерной (2Б) постановке задачи (система уравнений «мелкой воды»). Считается, что результаты, полученные с помощью таких методов, вполне достоверны и надежны.
Это действительно так, но только применительно к простым случаям (большая протяженность области решения, отсутствие в нижнем бъефе препятствий и сложной топографии). Известно, что уравнения «мелкой воды» (в приближении Буссинеска, Сен-Венана) обладают недостатками, обусловленными упрощениями, с помощью которых они получены из уравнений Навье-Стокса, и не способны полностью описывать движение волны прорыва по пойме реки.
Задачи, связанные с определением
ной близости от прорана, в руслах со сложной геометрией, а также при взаимодействии волн с преградами и сооружениями на пойме, могут решаться только с помощью численного моделирования в трехмерной (по пространству) постановке. Использование математической модели, основанной на 3Б эволюционных уравнениях Навье-Стокса, с учетом интенсивно изменяющейся свободной поверхности в режиме турбулентного течения, предъявляет очень серьезные требования к математическим алгоритмам нахождения решения, методам аппроксимации и вычисления. Выполнение 3Б моделирования волн прорыва стало возможно только с недавнего времени, когда появились достаточно мощные вычислительные машины и алгоритмы, позволяющие решать трехмерные уравнения гидродинамики за вполне разумное время.
Для численных расчетов авторами был применен метод, основанный на решении трехмерных эволюционных уравнений Навье-Стокса, ранее представленный и описанный в работе [1], а также
алгоритм [2], адаптированный на параллельную среду, с использованием графического процессора компании NVIDIA технологии CUDA [http:/developer.download.nvidia.com/ compute/cuda/2_3/toolkit/docs/ NVIDIA_CUDA_Programming_Guide_2.3 .pdf].
Применение данной технологи позволяет повысить скорость расчета приблизительно в 180 раз, что, безусловно, делает реальным и целесообразным использование трехмерных по пространству уравнений при расчете параметров волны прорыва.
Результаты расчета волны прорыва по данному методу сопоставлены c данными физического эксперимента [3]. Для убедительности сопоставление результатов 3D моделирования выполнено в плоскости действия вектора гравитации.
Физический эксперимент проводили в бассейне с «мокрым» руслом в нижнем бьефе (дно - ровное, стеки - твердые, непроницаемые, вертикальные), разделенном в начальный момент времени непроницаемой перегородкой (затвор), создающей перепад уровней. Схема опытной установки представлена на рис. 1. В ходе экспериментов определялись параметры волнообразования и последующего течения, возникающего в результате мгновенного образования прорана.
Для численного эксперимента была создана компьютерная модель опытной установки, описанной в [3]. Геометрия области расчета построена по точным техническим характеристикам эксперимен-
Затвор
Рис. 1 Схема экспериментального гидравлического лотка. Ширина по оси у - 0,27 м
тальной установки. Все начальные условия при моделировании были приняты такие же, как и в лабораторных экспериментах.
При дискретизации области решения использовалась сетка тетраэдров, представляющих собой связанные объемные фигуры, покрывающие всю область решения (рис. 2). Для генерации сетки использовался Бе1аипау алгоритм, позволяющий создавать сетку, равномерно изменяющуюся по объемом фигур [4, 5].
Для моделирования турбулентного течения применена модель динамики больших вихрей, более детально описанная в [6]. Минимальный и максимальный размеры элементов сетки вычисляются таким образом, чтобы вмещать все масштабы движения до инвариантного участка переноса энергии турбулентного течения (рис. 3).
На рис. 3 приведены результаты измерений в различных течениях спектральной плотности энергии Ек пульсаций продольной составляющей скорости ук'.
Нижний бьеф
Место установки затвора
ж vi X
Рис. 2. Трехмерная модель экспериментального лотка, вид с верхнего бьефа. Область расчета покрыта неструктурированной адаптивной сеткой тетраэдров
№ 2' 2011
(39
106 105 104 103 ю2 10 1 10" 10" 10" 10" 10"
¡Чл х-1
------ /_ ! + -2 а-3
________ (л ----V5 к о-4 0-5
------ ----ооо о\ ч-6 4-7
,, ___ 7 ^^тУ
№ 1
!
0 \
10"5 10"
10"
10"
10"
1
Рис. 3. Спектры плотности энергии пульсаций продольной составляющей скорости для различных течений: 1 - приливно-отливной канал; 2 - круглая струя; 3 - течение в трубе; 4 - течение с продольным сдвигом; 5 - след за цилиндром; 6 - турбулентность за сеткой; 7 - пограничный слой
По оси ординат отложена безразмерная спектральная плотность Ек, определяемая соотношением
по оси абсцисс - безразмерное волновое число к п,
где к1 - волновое число, определяемое как параметр, обратный линейному масштабу струк-
Г з Л1/4 V
туры (вихря); Л —
- микромасштаб тур-
булентного течения, по А. Н. Колмогорову; V -коэффициент кинематической вязкости жидкости; е - скорость диссипации турбулентной энергии.
На рис. 3 видно, что структура крупномасштабных вихрей (малые меняется в зависимости от типа течения и числа Рейнольдса, а мелкомасштабные вихри, участвующие в диссипации (большие к1п) являются универсальными, не зависящими от типа течения и числа Рейнольдса. Крас-
ной вертикальной линией отделен инвариантный участок масштабов волновых чисел от масштаба, зависящего от макроскопических параметров геометрии и кинематики течения. Следовательно, усреднение на масштабах к^ц = 2,0-Ю1 является инвариантным относительно групп переноса и поворота, что делает такое усреднение независящим от геометрии и кинематики любой моделируемой задачи.
Этот факт позволяет выбрать простейшую модель замыкания - так называемую подсеточную модель. Такой подход полностью исключает зависимость выбранной подсеточной модели от каких-либо индивидуальных свойств моделируемой задачи.
Размер разрешения пространственных и временных интервалов для дискретной модели уравнений Навье-Стокса должен быть строго определен. Очевидно, что этот масштаб зависит от числа Рей-нольдса, определенного как
д_Ш _4ЁНН ^
V V '
где и = ^Н - характерная скорость течения
жидкости; Ь = Н - характерный масштаб движения и выбранного масштаба усреднения по волновым числам.
Известно, например, что скорость диссипации энергии для полностью развитого турбулентного режима (Ие > 20 000) может быть определена так [7]:
тъ
ил
8 =
тогда микромасштаб А. Н. Колмогорова для развитого турбулентного течения
Г| =
ч1/4
Яъиъ У
Учитывая, что волновое число связано с масштабом движения как к = 2п/п, для данной задачи получаем:
Г т Vм
ЛьЕБ _
2пх\ к{\1
2 п кр1
л3 и3
^ у
т. е. длина грани тетраэдра г|ьвд = 0,022511881 « 0,023 м.
Определенный масштаб линейного элемента пТЕЯ позволяет узнать число элементов, на которое необходимо разбить расчетную область. Так, для прямоугольной области с размерами Н = 0,2 м, В =
0,27 м, L = 10 м и объемом V = h-b-l = 0,54 м3 можно определить минимальное число элементов, зная объем одного элемента разбиения области. Поскольку используется сетка тетраэдров, то необходимо определить объем одного элемента:
Ft^tiLB =1,43389-10«.
Таким образом, число элементов со-V
ставляет п = — = 376 598 шт. FT
Далее к Уравнениям применяют ту или иную модель подсеточного замыкания масштабов усреднения. В данной работе используется подсеточная динамическая модель Лэре, которую применяли во многих работах по математическому моделированию турбулентных течений несжимаемой жидкости [6].
Во время численных экспериментов был изучен процесс формирования волны прорыва в начальной стадии движения. Обработка полученных результатов показала следующее: после разрушения плотины массы воды верхнего бьефа
образуют волну опрокидывания, из-за чего возникает аэрированная пробка воды, которая, находясь в теле волны, вместе с ней продвигается далее по лотку.
Проверка адекватности математической модели выполнена сопоставлением результатов численного и физического экспериментов [3]. На рис. 4 изображены свободные поверхности воды при распространении волны в одинаковые моменты времени.
Обнаружено совпадение полученных в результате численного и физического экспериментов внутренней структуры тела волны, формы свободной поверхности и пройденного фронтом волны пути в моменты времени 0,156; 0,219; 0,281; 0,343; 0,406; 0,468 и 0,531 с. Такое совпадение свидетельствует о высокой точности применяемого метода для расчета параметров волны прорыва на расстоянии до 10 напоров.
Авторам не удалось найти опубликованных результатов натурных или лабораторных исследований волн прорыва, доступных для сопоставления и рассмотрения диапазонов адаптации применяемого
0,156 с
<Ъ|
0,219 с
0,281 с
0,343 с
0,406 с
0,468 с
I ■ Ч
0,531 с
_
0,156 с
0,219 с
0,281 с
с
0,343 с
0,406 с
С^ч
0,468 с
0,531 с ШГ С___
а б
Рис. 4. Сопоставление результатов физического (а) и численного экспериментов (б). Продольное сечение, после затвора. Изоповерхности свободной поверхности волны в различные моменты времени
№ 2' 2011
(41)
1
2
3
4
5
6
7
метода численного 3Б моделирования.
Выводы
Анализируя полученные результаты, можно утверждать, что использовать численные методы, основанные на 3Б уравнениях Навье-Стокса, можно при исследованиях распространения волн прорыва в реальных топографических условиях.
1. М. А. Волынов, И. В. Гугушвили, Н. М. Евстигнеев. Применение численных методов интегрирования трехмерных нестационарных уравнений гидродинамики при расчете распространения волны прорыва // Природообустройство. -
2009. - № 5. - С. 75-80.
2. Евстигнеев Н. М. Численный метод решения уравнений Навье-Стокса на неструктурированных сетках с применением Лагранжево-Эйлерового метода // Научно-технические ведомости СПбГПУ. -
2010. - № 1 (93). - C. 163-170.
3. Janosi I. M., Jan D., Szabo K. G. and Tel T. Turbulent drag reduction in dam-break flows. Experiments in Fluids 37. - 2004. - P. 219-229.
4. Cignoniz P., Montaniz C., Scopigno R. DeWall: A Fast Divide & Conquer Delaunay Triangulation Algorithm in Ed //
The Computer J. - 2006. - № 19(2). -P. 178-181.
5. Su P. and Scot Drysdale R. L. A
comparison of sequential delaunay triangulation algorithms // 11th ACM Computational Geometry Conf. Proc. -Vancouver, Canada: ACM Press, 1995. -P. 61-70.
7. Chunlei Liang, Evstigneev N. A study of kinetic energy conserving scheme using finite volume collocated grid for LES of a channel flow: Proceedings of the international conference on numerical methods in fluid dynamics. - King's College London, Strand, WC2R 2LS. -2006. - P. 61-79.
7. Pearson B. R., Krogstad P. A. and Van de Water W. Measurements of the turbulent energy dissipation rate // Phys. Fluids 14. - 2002. - № 1288.
Материал поступил в редакцию 01.02.10. Волынов Михаил Анатольевич, кандидат технических наук, доцент, заместитель директора
Тел. 8 (499) 153-21-33
Гугушвили Иракли Викторович, младший научный сотрудник
Тел. 8 (499) 153-21-33, 8-926-015-05-50
E-mail: Gugushvili_I@mail.ru
