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

УДК 532.526, 532.529

Д.И. Попов, А.М. Сагалаков

Исследование устойчивости течений в канале цилиндрической геометрии спектральными методами

D.I. Popov, A.M. Sagalakov

The Spectral Method of Flow Stability Investigation in Cylindrical Geometry

Представлен псевдоспектральный алгоритм решения линейных уравнений в частных производных с использованием полиномов Чебышева. Проанализировано распределение энергии пульсаций монодисперс-ной смеси. Установлено, что спектр для смеси может состоять из непрерывной и дискретной частей.

Ключевые слова: Метод Петрова-Галеркина, псев-доспектральный алгоритм, параллельные течения, монодисперсная смесь, показатели Ляпунова.

The pseudo-spectral algorithm to solve linear PDE using Chebyshev polynomials is produced in this work. The authors consider the distribution of pulsation energy for the mono-dispersed mixture. It is found that spectrum can consist of continuous and discrete sets.

Key words: Petrov-Galerkin method, pseudo-spectral algorithm, parallel flow, mono-dispersed mixture, Lyapunov exponents.

Введение. Наше внимание будет приковано к спектральным схемам аппроксимации некоторого комплексного гильбертова пространства и выбора базисов, соответственных для особенности геометрии течения и условия несжимаемости. Последнее, как известно, накладывает определенные требования на базисные элементы и элементы двойственного базиса, попытка удовлетворить которым в трехмерном случае обусловливает дополнительные технические сложности (особенно в случае спектральных методов в «чистом» виде). Здесь базисные элементы, на основе которых строится приближение, не обязательно будут совпадать (с точностью до операции взятия комплексного сопряжения) с элементами двойственного базиса. В рамках схемы, больше известной как соленоидальный метод Пет-рова-Галеркина [1, 2], будут даны примеры соответствующих базисов.

Уравнения, описывающие устойчивость параллельных течений в трубе монодисперсной жидкости в методе нормальных мод, могут быть записаны следующим образом:

Щ+( XV) U + (UV) V =

= _Vp+ду / Я + (//т)(У,-VI),

^ +(V2V)U + (^), =(У^2 )/т,, (1)

^ =0, дf/ дt+1т^(f^)}=0,

VI =у, VI =у.

1 1дО ь 2п 1дО 25

Здесь величина т = Ж - безразмерное время скоростной релаксации; Я - число Рейнольдса; £ - параметр, определяющий степень дисперсности примеси; О - ограниченная область Е3 с кусочно-гладкой границей дО.

Рассмотрены профили сдвигового течения и (г) , соответствующие напорному течению в зазоре между коаксиальными цилиндрами и (г )=А • г2 - В • 1п г +С, где коэффициенты определяются из граничных условий, и течению Гагена-Пуазейля и (г)=1 - г2. Решение спектральной задачи рассматривается на векторах из L2(О) с компонентами класса W012(О)nW22(О)

(в случае V) и W21 (О) (для V). Для решения применялся метод Петрова-Галеркина.

Конструирование соленоидального базиса. Известно [3, 4], что пространство L2(О) может быть представлено прямой суммой подпространств, одно из которых является замыканием бесконечно гладких соленоидальных векторных полей с компактным носителем, а другое - градиентов однозначных функций. Для всюду плотного множества векторов из L2 (О) можно записать VV=a1VU1 +a2VU 2 +6^^ (далее суммирование

производится по немому индексу). Представим

иГ 1 к 2 к 3 к“|Т ^ ^

1 =1 а1ке1,а1ке2,а1кез I и воспользуемся линейной

независимостью базисных элементов, тогда условие несжимаемости даст следующее соотношение:

а,/ \)+аг / 2(-)+а3 / 3()=0,

где ГО+а,*^*(0+а^зез () .

Разрешая полученное уравнение относительно одного из неизвестных коэффициентов для каждого к, получим набор соленоидальных базисных векторов. В случае малых гармонических по однородным координатам возмущений условие соленои-дальности выглядит наиболее просто

а1(гкг)'/г +а2тм>ф /г+а3тм>ф=0 ,

(гКг)' / г-/т*ф / г-/ай^=0 , где а е (0, +ГО); т = 0,1,2,_- осевое и азимуталь-

ное волновые числа.

В работах [1, 2] приведен пример такого базиса, при этом задача решается псевдоспектральным методом. Несмотря на то, что подобные схемы дают заполненные матрицы с довольно высокими числами обусловленности, при сравнительно небольшом разрешении аппроксимации точность схемы сравнима с результатами, рассчитанными спектральными методами [5].

В случае напорного течения в кольцевом зазоре традиционно для конструирования пробного и поверочного базисов предлагается применять элементы кк =(1 -2)Т и ик =(1-?2)2Тк, где t е [-1,1], Тк - полином Чебышева. Выпишем соответствующие производные

К =(Н2) п-2 Т,

К =(Н2) Т"-4 Т-щ, и* =(1-2)-2Ы* ,

ик =(1- 2) К -4 Ч-2 К.

Ясно, что выбор базиса не является однозначным (например, можно использовать функцию вида ик =Т* -2(к+2) / (к+3)Тк+2 +(к+1)/(к+3)Тк+4). В настоящей работе мы руководствовались сравнительно простой реализацией вычислений при достаточной точности расчетов, поэтому выбор базиса был наиболее простым для тФ0

%к =а^1+а^^2 , у =[/тик;-(гик)';0]Т у2 =[0; -1агкк; /ткк ]Т, р =[-/тик; -(гик)';0]Т , <р2 =[0; /аг*к; -ткк ]Т ,

V2k =а1^1+а2^2 +а3^3 =ак Кк ^ .

Полиномы Чебышева, как и все ортогональные многочлены, могут быть рассчитаны с использованием трехэлементных соотношений, но для расчетов полиномов Чебышева можно применять и тригонометрические выражения, что может упростить алгоритм вычислений. Однако оказывается, что рекуррентные соотношения более надежны и устойчивы к погрешности округления, обусловленной разрядностью соответствующей арифметики. В данной работе для решения обобщенной задачи на собственные значения использовались библиотечные процедуры ЬЛРЛСК95. Для улучшения обусловленности результирующих матриц применялась вариационная формулировка задачи, позволяющая понизить порядок производных.

Результат действия оператора дифференцирования второго порядка на пробные функции может быть записан следующим образом:

(Д у1)1 =Д у1-уI / г 2 -(2/т / г 2 у2,

(Ду2^ =-(2/т / г ^у^,

(Д у)2 =Д у2 -у2 / г 2 +(2/т/г 2 у1, (2)

(Д у2 )2 = Д у22-у22 / г 2,

(Ду)з =0, (Д у2 )з =Ду^.

где Д^=(1/г)(г^')'-(т2/г2)^-а2£ .

Для течения Гагена-Пуазейля характерным является то обстоятельство, что уравнения имеют геометрическую сингулярность на оси течения, что должно накладывать определенные условия на локальное поведение решения. Вообще говоря, задачи, когда уравнения имеют некоторую особенность, часто встречаются в практике. Общий способ регуляризации решения в окрестности сингулярности состоит в том, чтобы наложить определенное условие на локальное поведение решения - бихевиори-стическое условие [5]. Возможны различные способы регуляризации (особенно численной) решения в окрестности сингулярности. Например, в работе [6] предложено преобразование координат, позволяющее обеспечить гладкость решения на оси и размещать узлы квадратур или сетки в методе конечных разностей непосредственно в точке сингулярности. В работе [7] предложена спектральная схема аппроксимации для цилиндрической геометрии с использованием полиномов Чебышева и Лежандра. Для разрешения сингулярности при г = 0 ставится однородное условие, которое рассматривается в качестве главного краевого условия для случая тф0 . В случае т = 0 элементы пробного и поверочного базисов должны быть конечными в нуле. В случае соленоидальных векторных полей для регуляризации задачи в полюсе необходима аналитичность решения, что подробно обсуждается в работах [1, 2, 5, 8]. Для цилиндрической геометрии условие аналитичности при различных т дает теорема Приймака-Миядзаки.

Теорема (Приймак-Миядзаки): Для аналитического при г<е для некоторого е>0 векторного

поля ^г,0)=е/т0[уг(г^ + Ув(г')eв+vz(г] радиальная, азимутальная и аксиальная компоненты должны удовлетворять следующим условиям: уг =гкЕ (г), ve=гgE (г) (т=0); уг =гт1ИЕ (г),

у0= г1”1-1 £е (г) (т ^ 0); =г^ёЕ (г) (V т), где

^Е, £Е, ЬЕ - четные аналитические функции.

В работе [7] при разработке схемы спектральной аппроксимации производится преобразование координат вида г = (/ + 1)/2, где / е I = (-1,1), что вполне обосновано, когда спектральные коэффициенты рассчитываются аналитически. Однако в псевдо-спектральных схемах подобное преобразование может привести к необоснованной концентрации квадратурных узлов в окрестности нуля, что может ухудшить точность расчета скалярного произведе-

ния (например, в случае, когда решение заметно меняется в пристенной области [5]). Поэтому считается уместным продолжать решение в область г е (-1,0) таким образом, чтобы выражение, стоящее внутри скалярного произведения (•, •), было четным. В работах [1, 2] обсуждается, каким образом четность или нечетность решений можно использовать для упрощения спектральных матриц и сокращения интервала интегрирования г е I = (-1,1).

В настоящей статье спектральная схема аппроксимации рассматривалась на основе элементов, предложенных в [7]. Пробные и поверочные базисы в псевдоспектральной схеме должны были удовлетворять следующим условиям: аналитичность в окрестности нуля для несущей компоненты смеси, ограниченность - для взвеси; четность выражения, стоящего внутри скалярного произведения; элементы пробного и поверочного базисов для взвеси на стенке удовлетворяли однородному условию, которое считалось главным краевым условием. Таким образом, элементы пробного и поверочного базисов можно записать в следующем виде: у1=\_/тик;-(гик) ' ; 0]Т, у2=[0; -/аг2к* ; /тгк* ]Т,

рх=\_-/тгик;-(г2ик) ' ;0]Т, р2=[0;/аг3Кк;-/тг2Кк]Т для т = 1;

у =[/тг \; -(г \) ' ; 0]Т, у2 =[0; -/аг2 Кк; /тгмк ]Т,

р =[-/тгик; -(г X) ' ; 0] Т, р2 =[0; /агкк; -/ткк ]Т для нечетных т;

у =\_/тгик;-(г2и*) ' ;0]Т, у2 =[0; -/агъкк; /тг2Кк]Т,

р =[-/тик;-(гик) ' ;0]Т, Р2 =[0; /аг2Кк; -/тгмк ]Т для четных т.

Здесь к=2к', где к' е N .

Элементы базисов для второй компоненты смеси соответственно запишутся в виде

#1 =[гик ;0;0]Т, # = [0; гк* ;0]Т, # = [0;0;к*]т

г/1 = [и*;0;0]Т, % =[0;К;0]Т, % = [0;0; гк*]Т для четных т;

£ =[И*;0;0Г, £ = [0; wk ;0]T, £ = [0;0;rwJT,

Г = [ruk;0;0]T, ri2 =[0;rwk;0]T, r = [0;0;wk]T для нечетных m.

Для удобства применения спектрального метода необходимо первую производную радиальной составляющей поля скоростей для несущей среды по r выражать через базисные элементы, что позволит понизить старшую производную в операторе и несколько улучшить число обусловленности соответствующей спектральной матрицы. Приведенные базисные элементы можно использовать и в «чистой» спектральной формулировке. В этом случае достаточно громоздкие преобразования с использованием рекуррентных соотношений или тригонометрического представления полиномов Чебышева, дадут необходимые выражения для коэффициентов спектральной матрицы. Однако в случае описанного базиса подобные усилия не будут оправданы в смысле сочетания простоты элементов базиса и хорошей обусловленности результирующих матриц. В таблице приведено сравнение величины первой пристенной моды при а=1, m = 1, R = 9600 для разрешения аппроксимации N = 50 (результаты приводятся в соответствии с [1]).

Спектр. Известно, что спектр в случае монодис-персной смеси представлен точечной и непрерывной частями. Точечный спектр состоит из дискретного набора собственных значений и хорошо изучен для уравнений Навье-Стокса [3, 4]. Непрерывная часть спектра обусловлена уравнениями для примеси задачи (1) на подпространстве L2(Q), состоящем из градиентов однозначных функций. Можно показать, что в этом случае спектральная задача может быть приведена к соответствующей задаче для оператора умножения на функцию ограниченного изменения, спектр которого непрерывен. Раздвоение «хвостика» обусловлено недостаточностью разрешения аппроксимации, причем данный эффект не устраняется повышением разрешения N, причем это наблюдается в точках численного спектра, расположенных практически на одной прямой Y = const, к которой асимптотически приближается точечный спектр.

Первая пристенная мода при a=1, m = 1, R = 9600

Л1 = -0.01317 + 0.95048 і (Salwen et. al., 1980)

Л1 = -0.013170795764 + 0.950481396668 і (Leonard & Wray, 1981)

Л1 = -0.0131707957650041151055 + 0.9504813966699031794843 і (Priymak & MLyazakL, 1998)

Л1 = -0.013170795764 + 0.950481396670 і (Meseguer & Trefethen, 1000)

Л1 = -0.0131707957 + 0.9504813966 і N = 50

На рисунках 1 и 2 представлена часть спектра для /= 0.1, а=1, Я=3 -103 при т = 1 и т = 2 (непрерывный спектр не отображен). Из рисунков видно, что изменение дискретных собственных чисел вполне согласуется с представлениями о трансформации спектра при относительно ограниченном

возмущении оператора [9]. В качестве невозмущенного оператора принимается оператор, соответствующий уравнениям Навье-Стокса. При этом рисунки иллюстрируют, что возмущение определяется не только параметрами примеси, но и азимутальными волновыми числами.

0.8 X 1

Рис. 1. Спектр в зависимости от степени дисперсности для f = 0.1, а=1, Я=3 • 103 при т = 1

0

-0.2

-0.4 У -0.6

-0.8

-1

0.2 0.4 0.6 0.8 X 1

Рис. 2. Спектр в зависимости от степени дисперсности для f = 0.1, а=1, Я=3 • 103 при т = 2

■гн ^ (□ о -0%- ■ ° * О 8 х о | ! ° % ! % ! О ^ п о» 08 О» о й ----“-«‘•--ш

1 О ь \ ? « |

□ -с1еаг . х -Ю“5 О -5-10-4 1

■ 1 1—■—1- — 1 ■ 1

с ё

Рис. 3. Распределение энергии пульсаций по сечению канала приf = 0.1, %=1 и £=2 -10 5 . 0 - суммарная энергия ч.ж.; 1 - г1(у) ; 2 - т2 (у) ; 3 - Р(у); 4 - (у) ;

5 - суммарная энергия смеси. а, Ь - т = 2; с, й - т = 3 Я=10 . (для ч.ж. а=1.5 , Я=104)

Ер

0.5'

-0.5

А

То і і г j і Т \ 1 J с 1 1 , 1

~ЬГ" чЖ/ J W Ті \ Т2 / ✓ Г

' • ■ —— т

Ер

0.5

-0.5

Ті , 1 * і /74 N f / \ / \ / \ / \ / \ / \ \ \ \ \ \ \ Vo

/ / / / / / / / / / . / t / * / / / / / / \ \ \ \ \ \ \ \ \ \ \ \ \ \ / / 1

'\DS \Гг A ' \ 4 4 Q ' Л '/ / / / / /

и A \ -s ■ 4 \ V > АҐ / V, ' / / r t / Do / / / • ■

0.5

0.6

0.7

0.8

0.9

0.1

0.2

0.3

0.4

Рис. 4. Распределение энергии в сечении канала для/ = 0.1, а=1, Я=104, т = 1, £=10 4, С = 0.260994069280 - 10.0633086565, нолик относится к чистой жидкости Сс1 = 0.2737887094 - 10.0472321996; а - пристенная мода; Ь - приосевая

силы Стокса, Б=г1+В - локальный избыток гене-

Распределение энергии в сечении канала.

Уравнения, описывающие баланс энергии пульсаций, несложно получить в общем случае, умножая уравнения (1) на сопряженные к полю скоростей и интегрируя по О . Здесь примем такие обозначения: гк - характеризуют обмен энергией между основным профилем и возмущениями для компонентов смеси; Б - вязкая диссипация; +Б2 - работа

рации энергии над диссипацией для первой компоненты смеси. Для случая коаксиальной геометрии существенно разделение возмущений на пристенные и «приосевые». Это проиллюстрировано на рисунке 3.

На рисунке 4 представлено распределение в случае течения Гагена-Пуазейля.

Библиографический список

1. Meseguer A., Trefethen L.N. A spectral Petrov-Galerkin formulation for pipe flow I: Nonlinear transitional stages // Oxford University Computing Laboratory, Numerical analysis group, technical report no. 00/18, september, 2000.

2. Meseguer A., Trefethen L.N. Linearized pipe flow to Reinolds number 1.E+7 // J. Comp. Physics. - 2003. -Vol. 186.

3. Темам Р. Уравнения Навье-Стокса. Теория и численный анализ. - М., 1981.

4. Ladyzhenskaya O.A. The mathematical theory of viscous incompressible flow. - New York, 1963.

5. Boyd J.P. Chebyshev and Fourier spectral methods. DOVER Publications, Inc., 2000.

6. Kamran Mohseni and Tim Colonius Numerical Treatment of Polar Coordinate Singularities // J. Comp. Physics. - 2000. - Vol. 157.

7. Shen J. Efficient spectral-Galerkin methods IV. Spherical geometries // SIAM J. Sci. Comput. - 1999. - Vol. 20, No. 4.

8. Canuto C., Hussaini M.Y., Quarterni A. & Zang T.A. Spectral methods in fluid dynamics // Springer series in computational physics. - Berlin; New York, 1988.

9. Като Т. Теория возмущения линейных операторов. - М., 1972.