CC BY
38
Вычислительные технологии
Том 18, № 3, 2013
Численное моделирование устойчивости локализованных возмущений в течении Пуазейля*
А. В. Проскурин, А.М. САГАЛАКОВ Алтайский государственный университет, Барнаул, Россия
e-mail: k210@list.ru
Предложен численный метод исследования устойчивости течений по отношению к локализованным возмущениям. Метод основан на решении задачи на собственные значения для линейных уравнений в частных производных. Аппарат Я-функций позволяет строить приближённые решения для областей почти произвольной геометрии, точно удовлетворяющие граничным условиям. Алгебраическая задача на собственные значения решалась итерационными методами.
Ключевые слова: гидродинамическая устойчивость, течение Пуазейля, уравнения с частными производными.
Введение
Исследование устойчивости и бифуркаций течений представляет собой сложную проблему. Распространённые методы прямого численного моделирования течений дают результаты, близкие к экспериментальным. Такие вычисления сложны, а их результаты не всегда надёжны. Анализ устойчивости течений дополняет методы прямого численного моделирования и часто бывает полезен для изучения переходных процессов в течениях.
Классическая постановка задачи устойчивости течений вязкой жидкости приводит к задаче Орра — Зоммерфельда [1] — краевой для системы обыкновенных дифференциальных уравнений. Это возможно благодаря наличию однородных направлений, относительно которых возмущение можно представить в виде элементарных волновых решений с определённым набором волновых чисел, например, для течения в плоском канале, в круглой трубе, между коаксиальными цилиндрами. С другой стороны, есть много течений, в которых картина развития возмущений вдоль или поперек потока неоднородна: течение в прямоугольной трубе, обтекание крыла и цилиндра и др. В таких случаях задача устойчивости сводится к краевой задаче для уравнений с частными производными.
Многолетние исследования показали [1, 2], что эффективное решение задач устойчивости течений вязкой жидкости возможно только численно. Для решения задач устойчивости, как правило, необходимо использовать численные методы высокого порядка точности. В одномерном случае можно также применять методы низкого порядка, достигая высокой точности за счёт измельчения расчётной сетки. В двумерном случае это уже затруднительно. Например, в работе [3] показаны преимущества метода конечных элементов высокого порядка по сравнению с обычным методом конечных элементов.
* Работа выполнена при финансовой поддержке федеральной целевой программы "Научные и научно-педагогические кадры инновационной России", контракт 14.740.11.0355.
Решение краевой задачи для уравнений в частных производных существенно зависит от формы области и граничных условий, поэтому численный метод должен учитывать геометрическую информацию о задаче. В методе конечных элементов исследуемая область задается множеством меньших областей — конечных элементов. В работе [4] обсуждается принципиально иной подход. При помощи аппарата Я-функций строится уравнение границ исследуемой области. Это уравнение используется для построения некоторой структуры (пучка функций), которая приближённо представляет решение краевой задачи с учётом граничных условий. Неизвестные параметры данной структуры можно найти стандартными численными методами: методом коллокаций, Бубнова— Галеркина, наименьших квадратов.
Рассмотрим бесконечно малое возмущение течения Пуазейля, локализованное по длине канала. Если амплитуда возмущения нарастает в фиксированной точке, неустойчивость течения называют абсолютной (см. [5]). Возможны случаи, когда возмущение нарастает, смещаясь вниз по потоку, а его амплитуда в неподвижной точке убывает. Такую неустойчивость называют конвективной. Кроме того, могут образовываться "турбулентные клубы" — локализованные в пространстве турбулентные структуры. В работе [6] при некоторых числах Рейнольдса в круглой трубе наблюдались перемежаемые течения, "в которых локализованные турбулентные структуры, окружённые практически ламинарными участками течения, сносятся вниз по потоку, сохраняя свои пространственные размеры". В связи с вышеизложенным представляет интерес исследование устойчивости локализованных возмущений в течении Пуазейля, что и является целью настоящей работы.
1. Постановка задачи
Рассмотрим плоское течение Пуазейля: течение вязкой жидкости между бесконечными параллельными плоскостями. Под действием постоянного градиента давления жидкость движется с постоянной скоростью. Направим ось х вдоль направления движения жидкости, а ось у — перпендикулярно плоскостям (рис. 1). Уравнение Навье — Стокса для функции тока имеет вид
д Дф _ д^ д+ д^ дД^ = —АЧ (1)
дЬ дх ду ду дх Яе '
Плоскость у
_Ь и = 1 _ у2 1 п 0 ь
-1 X
Плоскость
Рис. 1. Конфигурация потока
где ф — функция тока, Ле = V0d/v — число Рейнольдса, V — максимальная скорость течения, й — полуширина зазора между плоскостями, V — коэффициент вязкости. В соответствии с методом линеаризации [7] решение (1) находим в виде
ф = фо + ф(х,у)ес (2)
где ф0 = у — у3/3 — стационарное параболическое решение, ф(х,у)ес* — малое возмущение (см. также [3, 8]), ф(х,у) — амплитуда, С определяет декремент затухания возмущения. Подставляя (2) в (1), получим
С(фхх + фуу) Т> (фхххх + 2фххуу + фуууу) 2фХ (1 у )(фххх + фуух) , (3)
ле
где нижние индексы означают производные.
Зададим условия на границе расчётной области П (см. рис. 1). На неподвижных стенках возмущение скорости обращается в нуль, на входе в область П возмущение задано равным нулю, на выходе из П, следуя [9], примем нулевой градиент скорости вдоль оси х. В этом случае возмущение может покидать расчётную область с ненулевой амплитудой. В работах [3, 5, 8], наоборот, использовались нулевые граничные условия. Это позволило упростить вычисления. Возмущения в данном случае не могут покидать область расчётов и остаются локальными. В работе [5] сделан вывод о том, что расчётная область при использовании граничных условий второго типа должна быть длиннее. В настоящей работе использовались граничные условия
ф = фх = фу = 0, (4)
соответствующие обращению возмущений в нуль. Таким образом, рассматриваются возмущения, локализованные в пространстве. Считается общепринятым, что линейная теория описывает возмущения только на начальном очень коротком отрезке времени, за который возмущение не успевает покинуть расчётную область.
2. Численный метод
Решение (3), удовлетворяющее граничным условиям, будем искать в виде
ф = ш ^ аз Тг(х)Тз (у^ (5)
где ш(х, у) — функция, определяющая уравнение границы ш = 0, Т(х), Т(у) — многочлены Чебышёва первого рода, а^ — неизвестные коэффициенты.
Определив множество точек коллокации в нулях полиномов Чебышёва
х* = Ь сов(п(2г — 1)/(п + 1)), г = 1,... ,п + 1,
у = С08(п(2^' — 1)/(к + 1)), э = 1,...,к + 1, получим задачу на собственные значения
Av = С Bv, (6)
где V = {а00, а01,... , ап(к-1), апк}. Матрица А будет содержать (пк)2 элементов. Полагая п = к = 60 — типичному значению для решения задачи Орра — Зоммерфельда, можно
оценить размер матрицы при использовании чисел двойной точности: (646) • 604 ~ 99 мегабайт. Добавление ещё одной пространственной координаты в разложении (5) равносильно умножению на п2, при этом необходимый для матрицы объем памяти становится слишком большим для персональных компьютеров.
Чтобы написать уравнение границы прямоугольника П, рассмотрим пересечение полос XI = 1 _ у2 > 0, х2 = Ь2 _ х2 > 0. Тогда искомое уравнение даётся, например, Я-конъюнкцией (см. [4])
Х1 Ло Х2 = (xi + Х2 - у x\ + x2)(x1 + x22)2 , (7)
где m — натуральное число. Основные вычисления были проведены при m = 0. Подстановкой x1 и x2 в (7) получим
ш(х,у) = 1 - y2 + L2 - x2 - V(1 - У2)2 + (L2 - x2)2. (8)
Производные от ш находим при помощи системы компьютерной алгебры Maxima. Выбор функции ш неоднозначен. Для проверки численного метода задачу (6) можно решать с использованием разных функций ш. В нашем случае применялись уравнение (7) при m > 0 и выражение
ш^,у) = ^2 - x2)(1 - у2). (9)
Первоначально задача на собственные значения (6) решалась с использованием QZ-алгоритма из библиотеки LAPACK, который успешно применялся нами в одномерном случае [10]. Однако вычислительные затраты этого метода оказались очень велики, а данный алгоритм применительно к нашей задаче иногда приводил к ошибочным результатам, поэтому для расчётов были использованы методы Крылова [11, 12], реализованные библиотеками ARPACK и SLEPc. Обобщённая задача на собственные значения (6) при помощи функций из библиотеки LAPACK сводилась к задаче AB-1v = Cv. Вычислялась часть спектра матрицы AB-1 с наибольшей действительной частью. Численные эксперименты показали эффективность обратных итераций с нулевым действительным сдвигом.
3. Результаты вычислений
Рассмотрим течение с L = 1.0 и n = k = 50. В табл. 1 приведены действительные части собственных значений C с наибольшей действительной частью при числах Рейнольдса от 10 до 1000, рассчитанные c помощью библиотеки LAPACK, SLEPC с использованием арифметики двойной точности и с помощью библиотеки NektarH—+, которая реализует метод спектральных элементов. Расчётная сетка для NektarH—+ была подобрана так, что погрешность приведённых в табл. 1 значений оценивалась величиной порядка 10-4 при Re = 1000. Для Re =10 все три собственные значения совпадают, но с увеличением числа Рейнольдса число совпадающих цифр уменьшается. При этом собственные значения, вычисленные с помощью LAPACK и SLEPC, близки друг к другу и отклоняются от значений, вычисленных с помощью Nektar++. При Re = 700 и Re = 1000 LAPACK определяет собственные значения, резко выделяющиеся из ряда и, по-видимому, не имеющие отношения к спектру задачи.
В табл. 2 приведены собственные значения с наибольшей действительной частью в зависимости от количества точек коллокации п и к, рассчитанные при помощи БЬЕРс. С увеличением п и к соответствующие собственные значения должны образовывать последовательность, сходящуюся к точному собственному значению. На практике после достижения определённой точности на результат могут влиять другие ошибки, свойственные конкретной реализации алгоритма. Как следует из результатов, приведённых в табл. 2, при И,е = 1000 действительные части собственных чисел сходятся по крайней
Таблица 1. Зависимость действительной части С собственных значений от числа Рей-нольдса при п = 50, к = 50, Ь = 1.0
Ие С Библиотека Ие С Библиотека
10 -3.1786 ЬЛРЛСК 400 -1.5169 ЬЛРЛСК
-3.1786 БЬЕРс -1.5169 БЬЕРс
-3.1786 -1.5156
50 -2.1144 ЬЛРЛСК 500 -1.4495 ЬЛРЛСК
-2.1144 БЬЕРс -1.4492 БЬЕРс
-2.1144 -1.4429 Nektar++
100 -1.9808 ЬЛРЛСК 600 -1.4048 ЬЛРЛСК
-1.9808 БЬЕРс -1.4039 БЬЕРс
-1.9806 -1.3850
200 -1.7517 ЬЛРЛСК 700 -0.3952 ЬЛРЛСК
-1.7517 БЬЕРс -1.3802 БЬЕРс
-1.7515 -1.3372 Nektar++
300 -1.6121 ЬЛРЛСК 1000 2.78053 ЬЛРЛСК
-1.6121 БЬЕРс -1.2963 БЬЕРс
-1.6120 -1.2306 Nektar++
Таблица 2. Зависимость собственных значений с наибольшей действительной частью от п и к при Ие = 1000, Ие = 5000, Ь = 1.0. Строки {1}, {2} соответствуют вычислениям с помощью Л-конъюнкции (7) при т = 1 и (9)
Ие = 1000 Ие = 5000
п к Действительная Мнимая Действительная Мнимая
часть С часть С часть С часть С
50 50 _1.293264539871856 1.087462148640639 1.942701048 1.559147236
60 60 _1.255768395039908 0.8781136507171899 0.6587661583 2.010687635
80 60 _1.230050190635487 0.8971942939279334 _0.4974945838 2.022046912
100 60 _1.230568610422895 0.8984693511654026 _0.8625281770 0.8222072725
100 70 _1.230568610049606 0.8984693516468647 _0.8625283592 0.8222064175
120 60 _1.230610310568465 0.8984959713463586 _0.8264486429 0.5237272589
140 60 _ 1.230611498906538 0.8984963168534776 _0.8244197295 0.5401157092
160 60 _1.230611522477419 0.8984963191778784 _0.8263675375 0.5437300819
180 60 _1.230611522273824 0.8984963207949902 _0.8272114112 0.5444320117
220 60 _1.230611522084476 0.8984963194842599 _0.8275200538 0.5445807151
220 70 _1.230611522758129 0.8984963207076385 _0.8275200601 0.5445807080
400 60 _1.230611518844661 0.8984963165508606 _0.8275373812 0.5445849618
{1} 220 60 _1.230613806678481 0.8985043165368143 _0.8275117537 0.5446018155
{2} 220 60 _1.230611536609566 0.8984964473559163 _0.8275529056 0.5445946517
мере до седьмого знака, а при Ле = 5000 — до четвёртого знака после запятой. Для контроля приведены также собственные числа, полученные вычислением с использованием Я-конъюнкции (7) при ш =1 и (9) (см. две последние строки в таблице).
Из постановки задачи следует, что собственные функции можно разделить на монотонные (У = 0) и периодические (У = 0). Собственные значения (С = X + гУ), соответствующие этим двум типам, при Ле = 1000, Ь =1 представлены на рис. 2. На рис. 3 показаны линии уровня действительной и мнимой части собственной функции, соответствующей собственному значению С с наибольшей действительной частью, при И,е = 1000, Ь = 1. В данном случае наблюдаются симметричные вихри. Аналогичный график приведён на рис. 4 при И,е = 1000, Ь = 4 в виде трёхмерной поверхности. Это монотонно затухающий вихрь с осью вблизи оси канала.
На рис. 5 представлены зависимости наибольшей действительной части собственных значений от Ь для разных чисел Рейнольдса от 103 до 5 • 103. Видно, что для каждой
У 30 20 10 0 -10
2.5 -2 -1.5 X
Рис. 2. Зависимость У от X при И,е = 1000, Ь = 1
Рис. 3. Линии уровня действительной (а) и мнимой (б) части собственной функции при И,е = 1000, Ь = 1; собственное значение С = (—1.230611522273843, 0.8984963207949633)
Рис. 4. Действительная часть собственной функции при Ие = 1000, Ь = 4; собственное значение С = (-0.1786773045527003, 0)
Рис. 5. Зависимости наибольшей действительной части собственных значений от Ь при Ие = 103 (1), Ие = 3 ■ 103 (2), Ие = 5 ■ 103 (3)
ветви величина X возрастает с увеличением Ь. В свою очередь сами ветви расположены в определённом порядке: при увеличении числа Рейнольдса значения X возрастают.
Заключение
Результаты проведенных расчётов позволяют заключить, что локализованные возмущения устойчивы в рассмотренных диапазонах параметров И,е и Ь.
Одна из задач данной работы — разработка численного метода для решения задачи гидродинамической устойчивости (3). Как уже отмечалось, эта задача содержит
малый параметр 1/Re, что значительно усложняет её решение и требует высокого качества приближённого представления решений. Предложенный численный метод алгоритмически проще многих других методов и учитывает гладкость решения. При L = 4, Re = 5000, n = 280, k = 60 вычисления на процессоре Phenom 2.2 ГГц занимают около 40 мин и примерно 7 Гб оперативной памяти. Продемонстрирована работоспособность метода вплоть до Re = 5000. Для сравнения — в одномерном случае аналогичный метод позволяет работать при Re ~ 107 [10].
В качестве примера рассмотрена область предельно простой геометрии — прямоугольник, однако алгоритм позволяет исследовать области более сложной геометрии. При этом граничные условия удовлетворяются точно. Ограничения могут быть связаны с громоздким видом граничной функции ш, производных от неё и с соответствующим ухудшением свойств пучка (5).
Авторы благодарят чл.-корр. РАН В.В. Пухначёва за полезные замечания. Список литературы
[1] Henningson D.S., Schmid P.J. Stability and Transition in Shear Flows. New-York: SpringerVerlag, 2001.
[2] Гольдштик М.А., Штерн В.Н. Гидродинамическая устойчивость и турбулентность. Новосибирск: Наука, 1977. 368 с.
[3] Gonzalez L., Theofilis V., Sherwin S.J. High-order methods for the numerical solution of the BiGlobal linear stability eigenvalue problem in complex geometries // Intern. J. for Numerical Methods in Fluids. 2011. Vol. 65, iss. 8. P. 923-952.
[4] Кравченко В.Ф., Рвачев В.Л. Алгебра логики, атомарные функции и вейвлеты в физических приложениях. М.: Физматлит, 2006. 416 с.
[5] Blackburn H.M., Barkley D., Sherwin S.D. Convective instabitity and transient growth in flow over a backward-facing step //J. Fluid Mech. 2008. Vol. 603. P. 271-304.
[6] Приймак В.Г. Волны и пространственно локализованные структуры в турбулентных течениях вязкой жидкости. Результаты расчётов // Матем. моделирование. 2010. Т. 22, № 2. С. 3-28.
[7] Юдович В.И. Метод линеаризации в гидродинамической теории устойчивости. Ростов-на-Дону: Ростовский гос. ун-т, 1984.
[8] Theofilis V. Advances in global linear instability analysis of non-parallel and three-dimensional flows // Progress in Aerospace Sci. 2003. No. 39. P. 249-315.
[9] Barkley D. Confined Three-Dimensional Stability Analysis of the Cylinder Wake. 2004. (http://arXiv.org/abs/physics/0405153v2).
[10] Proskurin A.V., Sagalakov A.M. Stability of the Poiseuille flow in a longitudinal magnetic field // Techn. Phys. 2012. Vol. 57, iss. 5. P. 608-614.
[11] Saad J. Numerical Methods for Large Eigenvalue Problems. Manchester Univ. Press, 1992. 346 p.
[12] HernAndez V., RomAn J.E., TomAs A., VidAl V. Krylov-Schur Methods in SLEPC. SLEPc Technical Report STR-7. Avaliable at: www.grycap.upv.es/slepc/documentation/ reports/str7.pdf. 2007. 13 p.
Поступила в 'редакцию 5 марта 2012 г., с доработки — 26 февраля 2013 г.
