Численное моделирование трехмерной неустойчивости обтекания короткого цилиндра Текст научной статьи по специальности «Физика»

Механика

УДК 532.5.011

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТРЕХМЕРНОЙ НЕУСТОЙЧИВОСТИ ОБТЕКАНИЯ КОРОТКОГО ЦИЛИНДРА

А. И. Алексюк1

В. П. Шкадова'

В. Я. Шкадов^

Плоскопараллельное обтекание бесконечно длинного кругового цилиндра становится трехмерным начиная с чисел Рейнольдса Re « 190. Соответствующую моду неустойчивости называют модой А. При Re « 260 в результате вторичной трехмерной неустойчивости (мода В) в следе возникают вихревые структуры с меньшим поперечным масштабом. В работе рассматривается процесс перехода к трехмерности для короткого цилиндра, ограниченного плоскостями. Длина цилиндра выбирается так, чтобы исключить неустойчивые возмущения моды А. Получены две моды неустойчивости, которые являются аналогами мод А и В, модифицированными под влиянием ограничивающих боковых плоскостей. Численные решения задач трехмерного обтекания строятся на основе уравнений Навье^ Стокса.

Ключевые слова: вязкая жидкость, трехмерные течения, обтекание цилиндра, неустойчивость, мода А, мода В.

The two-dimensional flow around an infinitely long circular cylinder becomes three-dimensional starting with Reynolds numbers Re « 190. The corresponding instability mode is called mode A. When Re « 260, structures with smaller cross-scale are formed in the wake as a result of a secondary three-dimensional instability (mode B). In this work, we consider the transition to three-dimensionality for a short cylinder bounded by planes. The length of the cylinder is chosen to eliminate the unstable perturbations of mode A. There have been found two modes of instability, which are analogues of modes A and В but modified under the influence of the limiting end planes. Numerical solutions of problems of three-dimensional flow are based on the Navier-Stokes equations.

Key words: viscous fluid, three-dimensional flows, flow around a cylinder, instability, mode A, mode B.

Введение. Плоскопараллельное обтекание бесконечного цилиндра неустойчиво относительно трехмерных возмущений при числах Рейнольдса Re > Re^ (Re^ и 190). Для конечных цилиндров периодическое течение со срывающимися параллельно оси цилиндра вихрями при Re < Re^ удается получить путем организации специальных условий на концах цилиндра [1, 2]. Другой способ получения параллельного следа связан с обтеканием длинных цилиндров, у которых относительные длины I = ^ (L, а — длина и радиус цилиндра) больше 1 ООО [3].

Бифуркация при Re = Re^ имеет субкритический характер: число Струхаля St скачкообразно падает на 5%, и в окрестности критического числа Рейнольдса наблюдается гистерезис [4]. Эту моду неустойчивости принято называть модой А. Масштаб развивающихся, периодических вдоль оси цилиндра структур около 8а. При Re > Rев (Rев ~ 260) в потоке формируются вихревые структуры меньшего масштаба 2а), что является следствием вторичной неустойчивости плоскопараллельного течения к трехмерным возмущениям (мода В).

Неустойчивые моды А и В и переход к трехмерности активно изучались экспериментально и теоретически, что продиктовано, в частности, стремлением к пониманию процессов развития вихревого следа и перехода к турбулентности в следе за телами, а также степени влияния на течение внешних воздействий [5-9]. Подробный обзор исследований задачи обтекания круговых цилиндров дан в работе [10].

В [11] на основе линейной теории Флоке получены значения критических параметров: Re^ = 188,5±1, А а = 3,96±0,02 для моды A, R ев = 259±2, Хв = 0,822±0,007 для моды В (Л — длина волны,

Алексюк Андрей Игоревич — канд. физ.-мат. наук, асснст. каф. аэромеханики и газовой динамики мех.-мат. ф-та МГУ, e-mail: aleksyukQmech.math.msu.su.

2 Шкадова Валентина Петровна — канд. физ.-мат. наук, вед. науч. сотр. НИИ механики МГУ.

3 Шкадов Виктор Яковлевич — доктор физ.-мат. наук, проф. каф. аэромеханики и газовой динамики мех.-мат. ф-та МГУ, e-mail: shkadovQmech.math.msu.su.

отнесенная к диаметру цилиндра). Показано, что для каждохх) Но € (Кед, 300) существуют два интервала неустойчивости (А™1П,А™ах), (А™11, Адах). Таким образом, в момент зарождения вторичной неустойчивой моды уже реализуется неустойчивая мода А. В экспериментальных исследованиях переход от моды А к режиму, где доминирует мода В, сопровождается плавным переносом энергии от волн первой частоты к волнам второй частоты [4|. В переходном интервале чисел Но наблюдаются две частоты.

При численном моделировании трехмерной неустойчивости на основе уравнений Навье Стокса возникают трудности, связанные с конечностью расчетной области. Это вынуждает ставить искусственные условия на ограничивающих цилиндр (с торцов) плоскостях. Такая постановка отличается от задачи обтекания бееконечнохх) цилиндра и может повышать устойчивость течения или наоборот приводить к неустойчивости. В экспериментальных исследованиях [12] рассмотрены подобные течения. Наблюдались двумерные обтекания цилиндра в горизонтальных и вертикальных пленочных потоках при 45 < Но < 560. Толщина пленки рмулировалаеь количеством хюдаваемохх) раствора и составляла 0,1 4- 10 мкм, диаметр обтекаемых цилиндров изменялся в пределах 1,2 4- 11 мм. Полученные в [12] числа Я! соответствуют плоекопараллельному обтеканию.

В настоящей работе рассмотрено обтекание короткохх) цилиндра, ограниченнемч) плоскостями с условиями отсутствия касательных напряжений и потока через эти границы. Длина цилиндра выбиралась из промежутка между двумя интервалами неустойчивых длин волн, т.е. I € (А^ах, А™1П). Цель исследований изучение отличий между процессами перехода к трехмерности для короткохх) и бееконечнохх) цилиндров, а также проверка предположения о том, что при длинах цилиндра I < А™1П и 11ел < Но < 11ев двумерное течение является устойчивым.

Постановка задачи и основные уравнения. Решение строится в области О (рис. 1), которая в декартовой системе координат (х1,х2,хз) задается следующим образом:

О = {х = (ж1,ж2,ж3) : х\ +х1 > 1; |ж3| < 1]ХВК < ад < Хвых; \х2\ < УСтк},

где Хвх, Хвых, 1'б„к константы, определяющие размеры области О. Граница области 90 состоит из входной границы Гвх = {х : х € 90п{ж1 = Хвх}}, боковых границ Г у = {х : х € {|а?21 = ^бок}} и Гг = {х : х е 90 П {|жз| = I}}, выходной храницы Гвых = {х : х € 90 П {х\ = Хвых}} и храницы тела Гт = {х : х € 90 П {х'1 + жI = 1}}- Размерные координаты связаны с безразмерными следующими соотношениями: ж' = аж,-.

МЖн

й"Ми',.?< > '»«' ;'уо

да» ЯЁ^йятщдб

■ Лй к 1%«8«8ГА »«ж ¿в

'ЯНяМ

к а ? ВАГ шп У' 1К 1 йР № ад 'ё ЛИ«) а№> 1 ИЖ1

£ Э 2С П > ; , \ / /'* - \ >8 * т 'ж •и мгуи 3

Рис. 1. Расчетная область и тетраэдральная сотка в сечениях = 0. хо = 0 и ж'з = 0

Задача сводится к численному решению трехмерной начально-краевой задачи для уравнений Навье-Стокса относительно неизвестного Y(^к,t) = (р,щ,и2,из,Т)*:

{ Л0У>4 + АгУ,г = (КгзУ!3 - Рг),г]

ОТ

щ = 0, —— = 0 при х € Гт; ап

щ = 5ц, Р = Роо, Т = тоо при х € Гвх;

ди

— = 0, Р = Роо,

Т = Too при х € ГЕ

дп

и • п = 0, (TijUjei) х п = 0, р = Poo, Т = тоо при х € TY

дТ

(1)

(TijUjei) х п = 0, — = 0 при х € Г2;

и • п = 0,

^щ = 5ц, Р = Роо, Т = Тоо при t = О, X е Q.

Здесь р, u = (ui,u2,us)*, Т — безразмерные давление, вектор скорости и температура; 7, Моо — показатель адиабаты и число Маха в набегающем потоке (в расчетах 7 = 1,4; Mqo = 0,1); Роо = JU ;

Too = —; e¿ — базис в системе координат Xi, п — нормаль к границе; (■) t = ^¡т', (•) к = ттг-',

7(.7—ijMoo ' ' к

по повторяющимся индексам ведется суммирование; (•)* — транспонирование. Матрицы Ао, Ai, Kij и вектор Р, (i,j = 1,2,3) возникают вследствие перехода от консервативных переменных к переменным давление-скорость-температура.

Для вектора теплового потока q'. тензора вязких напряжений r¿ • и внутренней энергии ё полагалось

q' = -kVT', T[j = /л (u'itj + u'jti - I 5iju'kt^j , ё = cvT' +

const,

где к = const — коэффициент теплопроводности, /л = const — коэффициент вязкости, су — удельная теплоемкость при постоянном объеме.

Обезразмеривание проводилось по формулам (размерные величины отмечены штрихами)

-и-

Ui =

Re =

Uoo'

2pODUODa

1

fx

P = —

Poo '

Pr =

AtCp к '

p =

St =

p

P00U00 2a

UooT'

T =

cvT'

и2 и2 :

^00 ^00

Moo =

Uo

2 = 7Pc*

00

Poo

где poo, Pod, Uoo — плотность, давление, скорость среды в набегающем потоке; ср — удельная теплоемкость при постоянном давлении; т — период вихреобразования; Рг — число Прандтля; Соо — скорость звука в набегающем потоке; ¿ = 1,2, 3.

Уравнения (1) описывают течения вязкого совершенного газа. При их решении число Маха было фиксированно: Моо = 0,1.

пт

Численный метод. Пусть Q = (J Те — разбиение расчетной области на тетраэдры; V, W —

е=1

пространства кусочно-линейных пробных и весовых функций соответственно, построенные стандартным для метода конечных элементов образом с учетом граничных условий Дирихле (1). Тогда стабилизированный метод конечных элементов GLS (Galerkin Least-Squares) для задачи (1) формулируется следующим образом: найти Y € V, такие, что для любого W € W

I (W • AqYj + W • AiYti + W,i • KijYj) dtt-

f пт Г

/ (W • KijYjm) dT + J2 (L*'W) • r(LY) dn = °>

г e=1Te

где т стабилизирующая матрица; L = Aq-щ + — дифференциальный оператор.

Здесь в стандартную формулировку метода Галеркина вводится еще одно слагаемое (последнее), которое вносит стабилизирующий вклад в уравнения. Метод является консервативным, а точное решение исходной краевой задачи удовлетворяет уравнениям в слабой форме. Детали метода описаны в [13].

Дискретизация частной производной но времени осуществлялась но тета-ехеме при 9 = 0,5 (ж = ^ од—I" 0(At2)). Дискретный аналог исходной задачи представляет собой нелинейную систему алгебраических уравнений, которая решалась на каждом шаге но времени итерационным методом Ньютона в сочетании с методом обобщенных минимальных невязок. Численный метод имеет второй порядок аппроксимации но времени и но пространству.

Построена параллельная реализация вышеописанного метода для систем с распределенной памятью на основе технологии MPI (Message Passing Interface интерфейс передачи сообщений). Распределение вычислений между MPI-процессами осуществлялось на основе разделения расчетной сетки на части [14]. На каждом процессе проводились вычисления с соответствующей частью узлов сетки и содержалась необходимая для счета информация о соседних (фиктивных) узлах, находящихся на другом процессе. Расчеты проводились на суперкомпьютерах СКИФ МГУ "Чсбышсв" и "Ломоносов" [15].

Результаты. Численные расчеты получены для цилиндров длиной i = 1,6 на расчетной сетке, изображенной на рис. 1: количество узлов и тетраэдров равно соответственно 305 029 и 1656 677, A"BX = —15, А"вых = 35, Y = 18. Узлы сетки сгущены вблизи тела и в области ближнего следа (шаг сетки 0,04 0,1). Выбранная для расчетов длина цилиндра удовлетворяет соотношению I G (А§ах, А™1П) для всех рассматриваемых чисел Re. Интервалы неустойчивых длин волн (значения А™1П, А™ах, А™п, А§ах) получены по теории Флоке в [11] и хорошо согласуются с экспериментальными данными.

Двумерные расчеты [8, 9] дают непрерывную зависимость St.(Re) в рассматриваемом диапазоне чисел Re (рис. 2) и хорошо согласуются с реальными течениями до чисел Re = Re^. Следует отметить, что данные по двумерному обтеканию [8, 9] (рис. 2) получены на основе метода, схожего с применяемым в настоящей работе, путем использования более подробных сеток. Отличие наблюдаемых в эксперименте частот от частот плоскопараллельного обтекания при Re > Re^ связано с трехмерностью реального течения. Экспериментальная зависимость St.(Re) имеет два скачка при Re = Re^ и Re = Пев и в окрестностях этих точек неоднозначна: при Re^ наблюдается гистерезис, при Reб одновременное существование двух частот в потоке [4|.

В интервале чисел Рейнольдса Re^ < Re < Re в в эксперименте реализуется мода А с характерным масштабом поперечных вихревых ctdvktvd около 8а. Ожидалось, что для этого режима при

I < А™1П обтекание короткого цилиндра также должно быть плоскопараллельным с частотой вихревой дорожки, соответствующей двумерному расчету. При числе Re = 200 данное предположение не подтвердилось: частота St. лежит между частотами, соответствующими двумерному и трехмерному обтеканию длинного цилиндра (рис. 2), а в ноле течения на рис. 3, а отчетливо видны вихревые структуры, имеющие трехмерный характер. Назовем эту моду неустойчивости модой А.

Течение при Re = 160 также обладает подобной трехмерной структурой, и его частота St. меньше, чем при двумерном обтекании (рис. 2). Таким образом, при Re = 160 плоскопараллельное обтекание короткого цилиндра неустойчиво: в потоке реализуется трехмерная мода А. Отметим, что для бесконечного

Рис. 2. Зависимость числа Струхаля = у—^ от числа Рейнольдса Т1с = 2/3ос^/о°а: сплошная линия двумерные расчеты [11]: круги эксперимент [10]: квадраты двумерные расчеты [9]: треугольники трехмерные расчеты, полученные в настоящей работе при I = 1,6: пунктирные линии значения Ие^ и Ш^в, полученные в [11]

цилиндра такое течение устойчиво. Очевидно, что плоекопараллельные решения задачи обтекания бесконечно длинного цилиндра также являются решениями задачи в усеченной области, рассматриваемой в работе. Но, как показывают результаты при Но = 160, усеченное решение может оказаться неустой чивым.

Характерная длина волны вихревых структур моды В при Но > Кед около 2а, что согласуется с результатами численного моделирования моды В для обтекания короткого цилиндра при Но = 280 (рис. 3, б). Однако ограничивающие торцевые храни искажают течение частота Я! для этого режима меньше, чем в эксперименте для длинного цилиндра (рис. 2). Назовем эту моду неустойчивости модой В.

а

■2 0 2 4 6 8 10 -2 0 2 4 6 8 10

х^ х^

Рис. 3. Обтекание короткого цилиндра при I = 1,6: а поверхности = ±0,1 (слева) и с^з = ±0,5 (справа) для моды А при Т1с = 200: б поверхности = ±0,5 (слева) и с^з = ±0,7 (справа) для моды В

при Т1с = 280

Структуры течения, соответствующие модам А, В, имеют различную пространственно-временную симметрию аналогично модам А и В [11]. У моды А знак продольной составляющей завихрен-

ности чередуется при Х3 = const, у моды В он постоянен. Картины течения на рис. 3, а качественно похожи на поля течений для моды А на полупериоде в поперечном направлении (см. [11, рис. 9, 10; 10, рис. 15]). В частности, в момент отрыва основного вихря (шз) вблизи его центра формируется продольный вихрь (wi). Как и для длинного цилиндра, при Re = 280 возмущения основного потока (рис. 3, б) в виде ненулевой продольной составляющей завихренности (wi) концентрируются между основными поперечными вихрями (соз). Для всех мод (А, В, А, В) трехмерность начинает существенно проявляться только на некотором расстоянии от цилиндра (рис. 3).

Заключение. Результаты работы показали, что в задаче обтекания короткого цилиндра, ограниченного плоскостями (с условиями отсутствия касательных напряжений и потока через них), существуют две моды трехмерной неустойчивости А и В, подобные модам А и В в задаче обтекания бесконечного цилиндра. Несмотря на то что моды А и А, В и В имеют схожие структуры, вопросу о возможности обобщения результатов численного моделирования мод А и В (перехода к трехмерности) на ограниченных цилиндрах должно уделяться особое внимание. Так, используемое в работе ограничение потока внесло дополнительную неустойчивость А при Re < Re^, что вызвало более ранний (по числу Re) переход к трехмерности в рассматриваемой постановке. Механизм формирования неустойчивой моды А требует более детального изучения.

Следует отметить, что предположение о том, что при длинах цилиндра I < А™111 и Re л < Re < Rед двумерное течение является устойчивым, не подтвердилось (для ограниченных цилиндров возникает неустойчивая мода А).

Работа выполнена при финансовой поддержке РФФИ (проекты № 14-01-31106, 12-01-00405).

СПИСОК ЛИТЕРАТУРЫ

1. Hammache М., Gharib М. An experimental study of the parallel and oblique vortex shedding from circular cylinder //J. Fluid Mech. 1991. 232. 567-590.

2. Miller G.D., Williamson C.H.K. Control of three-dimensional phase dynamics in a cylinder wake // Exp. Fluids. 1994. 18, N 1-2. 26-35.

3. Norberg C. An experimental investigation of the flow around a circular cylinder: influence of aspect ratio // J. Fluid Mech. 1994. 258. 287-316.

4. Williamson C.H.K. The existence of two stages in the transition to three-dimensionality of a cylinder wake // Phys. Fluids. 1988. 31. 3165-3168

5. Navrose, Meena J., Mittal S. Three-dimensional flow past a rotating cylinder //J. Fluid Mech. 2015. 766. 28-53.

6. Persillon H., Braza M. Physical analysis of the transition to turbulence in the wake of a circular cylinder by three-dimensional Navier-Stokes simulation //J. Fluid Mech. 1998. 365. 23-88.

7. Алексюк A.M., Шкадова В.П., Шкадов В.Я. Гидродинамическая неустойчивость отрывного обтекания кругового цилиндра вязкой жидкостью // Вести. Моск. ун-та. Матем. Механ. 2010. № 5. 51-57.

8. Алексюк А.И., Шкадова В.П., Шкадов В.Я. Возникновение, развитие и затухание вихревой дорожки в следе за обтекаемым телом // Вести. Моск. ун-та. Матем. Механ. 2012. № 3. 24-32.

9. Алексюк А.И. Исследование отрывных обтеканий тел методом численного решения уравнений Навье-Стокса: Канд. дис. М., 2013.

10. Williamson C.H.K. Vortex dynamics in the cylinder wake // Ann. Rev. Fluid Mech. 1996. 28, N I. 477-539.

11. Barkley D., Henderson R.D. Three-dimensional Floquet stability analysis of the wake of a circular cylinder // J. Fluid Mech. 1996. 322. 215-241.

12. Wen C.-Y., Lin C.-Y. Two-dimensional vortex shedding of a circular cylinder // Phys. Fluids. 2001. 13, N 3. 557-560.

13. Hughes T.J.R., Franca L.P., Hulbert G.M. A new finite element formulation for computational fluid dynamics. VIII. The Galerkin/least-squares method for advective-diffusive equations // Comput. Methods Appl. Mech. Eng. 1989. 73, N 2. 173-189.

14. Karypis G., Kumar V. A fast and high quality multilevel scheme for partitioning irregular graphs // SIAM J. Sci. Comput. 1998. 20, N 1. 359-392.

15. Воеводин В.В., Жуматий С.А., Соболев С.И., Антонов А. С., Брызгалов П.А., Никитенко Д.А., Стефанов К. С., Воеводин В.В. Практика суперкомпьютера "Ломоносов" // Открытые системы. 2012. 7. 36-39.

Поступила в редакцию 13.02.2015