CC BY
46
УДК 532.5
Эффективный спектральный метод для исследования устойчивости дисперсных течений
Д.И. Попов, Р.М. Утемесов
Алтайский государственный университет (Барнаул, Россия)
Effective Spectral Projection Method for Stability Analysis of Disperse Flows
D.I. Popov, R.M. Utemesov
Altai State University (Barnaul, Russia)
Рассмотрена разработка эффективного спектрального алгоритма Галеркина для эллиптических задач. В качестве приложения описанного метода рассмотрена задача устойчивости параллельного дисперсного течения. Интерес к высокоэффективным спектральным алгоритмам обусловлен современными расширениями промышленных языков программирования (например, C++14), которые предоставляют средства для создания высокопродуктивного и надежного параллельного кода. Изучение движения сред с особыми свойствами является одной из актуальных задач вычислительной и математической гидродинамики, в частности, исследование свойств линейных операторов динамических систем. Для разработки алгоритма применяются линейные комбинации полиномов Лежандра. Показано, что спектральные коэффициенты основных дифференциальных операторов могут быть рассчитаны по явным алгебраическим формулам. Приводится сравнение точности алгоритма с другими методами. Проанализирован спектр устойчивости течения Куэтта — Пуазейля дисперсной смеси. Установлена форма достаточных условий устойчивости. Численно обнаружены области метастабиль-ности для дисперсных течений.
Ключевые слова: метод Галеркина, дисперсная смесь, спектральная аппроксимация, гидродинамическая устойчивость, течение Куэтта — Пуазейля.
БОТ 10.14258Дгуа8и(2016)1-08
This paper is concerned with development of an efficient spectral Galerkin algorithm for elliptic problems. The stability problem of parallel disperse flow is employed as an application of described method. The interest in highly effective spectral projection algorithms is raised with development of contemporary extensions of industrial programming languages (C++14, for example), which grants frameworks for constructing high performance and robust parallel programming code. The exploration of flows with special properties is the problem of numerical and mathematical hydrodynamics. The central focus of exploration is to describe properties of linear operators of corresponding dynamical systems. The development of algorithm is based on Legendre Galerkin approach. It is shown, that spectral coefficients of differential operators are explicitly expressed with algebraic formulas. The test of algorithm accuracy is performed. The stability spectrum of disperse Couette-Poiseuille flow is analyzed. The form of sufficiency conditions of stability is established. The metastability regions of disperse flow are observed numerically.
Key words: Galerkin method, disperse flow, spectral approximation, hydrodynamic stability, Couette-Poiseuille flow
Введение. Двухфазные течения широко распространены в природе и приложениях. Под дисперсной смесью обычно понимают смесь газа и взвешенных в нем твердых частиц, смесь газа и капель, жидкости и пузырьков, жидкости и твердых частиц [1]. Примерами промышленного применения дисперсных течений могут послужить псевдоожиженный слой, пневмотранспорт, теплообменник с гранулирован-
ным теплоносителем, нефтегазовая промышленность, транспортировка угля и руд и т. д. (см. [1]). Во многих промышленных процессах требуется строгий контроль над режимами дисперсного течения. Механизм для предвидения переходов или изменений режимов течений предоставляет теория устойчивости [2].
Эффективные спектральные алгоритмы решения краевых задач для интегро-дифференциальных
уравнений по-прежнему являются одними из центральных предметов исследований вычислительной математики. Например, в работе [3] развивается проекционный метод для моделирования ламинарно-турбулентного перехода, обсуждается оптимизация вычислений. В работе [4] метод Галеркина и полиномы Лежандра применяются для решения интегральных уравнений Фредгольма-Хаммерстайна, установлена скорость суперсходимости приближенного решения.
Исследование ламинарно-турбулентного перехода в двухфазных средах осложнено рядом обстоятельств, таких как наличие многих масштабов, характеризующих задачу, тепломассообмен на границе раздела фаз, влияние параметров примеси на коэффициенты переноса. В обзорной работе [5] обсуждается актуальное состояние исследований турбулентности многофазных дисперсных течений. Исчерпывающий перечень литературы отражает достигнутые успехи. В частности, авторы [5] среди перспективных направлений для фундаментальных исследований выделяют изучение механизмов модуляции турбулентности [6] частицами примеси. При численном моделировании крупномасштабных эффектов в турбулентном потоке в случае больших частиц вполне приемлем двухско-ростной (двухжидкостный) подход [5, 6]. Для этого необходимо знать структуру собственных функций соответствующего линейного оператора [2].
Известно, что для дисперсной смеси достаточные условия устойчивости состояния равновесия принимают сложную форму [7-10]. Например, в работе [7] была рассмотрена линейная устойчивость течения Куэтта в случае, когда объемная концентрация примеси не являлась малой величиной. Было обнаружено, что течение может быть неустойчивым при неоднородном распределении концентрации облака частиц в канале [7]. В работах [8-10] численно обнаружены эффекты существенного повышения порога устойчивости для конвективных течений [10] и в случае изотермического движения [8, 9] дисперсной смеси, а также области метастабильности при закритических режимах течения [8, 9].
В данной работе рассмотрен эффективный спектральный метод решения первой краевой задачи для системы дифференциальных уравнений, описывающих устойчивость движения смеси.
1. Постановка задачи. Монодисперсная смесь — это гидродинамическая модель двухфазной системы, состоящей из вязкой жидкости и взвешенных в ней твердых частиц. Движение смеси можно описывать двухскоростными уравнениями механики сплошных сред, когда примесь рассматривается как формальный эффективно невязкий газ, средние нормальные напряжения в котором обращаются в нуль, а несущая среда является несжимаемой. Межфазное взаимодействие определяется силой Стокса.
При исследовании устойчивости течения к малым возмущениям линеаризация производится в окрестности стационарного решения. Например, в плоскопараллельном канале основное решение — профиль напорного течения — можно принять в виде
и (х2)=и (х2)ех =[(1 - А)(1- х22)+Лх2 ] е^.
Здесь Л е[0;1] — безразмерная константа; е1, е2 — орты декартовой системы координат, ориентированные в продольном и поперечном направлениях к течению. На стенках канала задаются условия непроницаемости и отсутствия проскальзывания для вязкой жидкости. Предполагается, что стационарный профиль течения и совпадает в обеих компонентах смеси и А\\и = 0.
Величина т = Б • К.е — безразмерное число, описывающее стоксовскую скоростную релаксацию, К.е — число Рейнольдса (далее V=1/К.е), Б — критерий, определяющий степень дисперсности примеси, р — относительная массовая доля примеси (величина р < 1 считается постоянной и определяется как отношение массы частиц к массе жидкости в единице объема). При этом ~ , Б~а2, где а — диаметр частицы.
2. Уравнения двумерных возмущений. В работах [1, 11] показано, насколько важно знать спектральный портрет линейного оператора соответствующей динамической системы. Здесь мы ограничимся изучением свойств спектра в моногармоническом приближении [12].
Пусть х = х1, у = х2, ех = е1, еу = е2. Для параллельных течений возмущение поля скоростей можно представить в виде гармонического колебания |й(х, у), (х, у)} = {и( у), у )}ехр(га х), с периодом, равным 2п / а, где а — соответствующее волновое число. Ограничимся двумерными возмущениями, учитывая преобразование Сквайра [12]. Условие несжимаемости для несущей жидкости позволяет задать функцию тока ф(х, у) такую, что йх =дф/ду, йу = = —дф/дх. Уравнения для амплитуд колебаний могут быть записаны следующим образом:
ХАф = ¡сЮДф — ¡аи "ф + +vA2ф + (р/т )в — (р/т )Дф,
Хпх = ¡аишх — и'п + (1/т)( — пх),
\пу = ¡аUwy — (1/т )(аф + пу), (1)
в = п:' — iаwy, ф(±1)=ф '(±1)=0,
(±1) = Пу (±1) = 0.
Здесь Д = й 2/йу2 —а2. Спектр задачи (1) описывает устойчивость стационарного двумерного течения двухфазной смеси. Присутствие собственного
числа А с неотрицательной действительной частью свидетельствует о неустойчивости. Обобщенное решение задачи (1) принадлежат множеству функций класса (см. [11])
ф е н2Ц)п N¿(7), ^,^^ е н'(I),
(±1) = ^^ (±1) = 0, I = (-1;1).
Из уравнений (1) следуют (см. [9]) достаточные условия устойчивости, которые можно записать в виде вариационного неравенства
Щ II н ц) I N1 н а)
р и н2
- —— Н
2^ "н(I)хн(I)
< v (+/ s)i «- н i; i ))
<
(2)
для любых {и, н} таких, что
1н (I)
и V(I) + РIHI н (I )хн (I) =
-а21 \ф\ 1н (i) + р\И н (i )Хн (I)=1
где
F(N, а) = 1Щ "IГ + 2а21\фII2
\Г> / || Г ||н(I) II Г ||н(I)
Y = max U'.
-а41\ф\
н (I)
= { е М(Ф^):Щ(± 1) = Ф'( ± 1) = о}},
^ _{ е М(Ф^):^к(± 1) = 0}}. i
Здесь через (Pk,Ph )w _ JPkPhwdy обозначено
(3)
ска-
лярное произведение в пространстве ^ (I). Тройку (I, ш, Ф н) называют спектральной проекционной системой [13].
Приближенное решение будет представлено парой [пн, wN} в виде рядов
ф* ( у)=£ *ТФк ( у),
+ UyNey =Е (0) ФкЧ - афк e7), (4)
■ WyN e y =£(ök1)^k +а<к^к ey ).
В результате получим обобщенную задачу на собственные значения:
Z <00'
Z(10) Z(20)
= А
Z(01) Z(11) Z(21) (00)
Z(02) Z(12) Z(22) 1(01)
(0)
(1)
(2)
(02)
Q() Q(01) Q
Q(10) q(11) q(12) q(20) q(21) q(22)
(0)
(1)
(2)
(5)
Здесь Z(lj),Q(lj) - матрицы размерности NxN,
я
« = [я^я^,...,я^}1', 0 < I,} < 2.
В качестве экземпляра класса Ф* можно использовать полиномы Лежандра [1к (у)}*=0 еФ*, где у е I, ш(у) = 1. Для построения наборов V*, используем многочлены следующего вида:
А = к -Ьк+2, 0 < к < * + 2,
фк _ Lk - 2
2к + 5
2к + 3
(6) (7)
Условием (2) определяется правая граница числового образа линейного оператора задачи (1).
3. Схема численного решения. Для построения эффективного спектрального алгоритма решения задачи (1) будем использовать метод Галеркина. Для этого рассмотрим на множестве I класс ортогональных с весом ш(у) > 0, у е I, алгебраических многочленов с рациональными коэффициентами:
Ф* ={Рк <Рк,Р I = ^ }}
Рассмотрим подмножества V* и линейной оболочки М(Ф*) множества Ф* такие, что
2к + 7 к+2 2к + 7 к+4 0 < к < N.
Очевидно, что многочлены (6), (7) являются линейно независимыми. Формула (7) может быть получена в результате следующей цепочки рассуждений:
Фк _ ßк" h+4 + ßf' Lk+2 + ß(3) Lk,
Ф' _ ß?L'k+i + (( + ßk3) ) - ßk3) (2k + 3)Lk+1,
ßk2) + ßk3) =-ßk1, ф' _ ßk^ (2k + 7)Lk+з - ßf (2k + 3)Lk+1. Таким образом, получим
-(1) = 2k + 3 „(2) = 2 2k + 5 ß(3) = 1
ßk = 2k + 7, ßk = 2 2k + 7, ßk =1,
что и требовалось.
В качестве спектральной проекционной системы примем тройку (I, w( y) = 1,{Lk }N+4). Коэффициенты матриц Z(j), Q(j) в уравнении (5) можно записать следующим образом:
^^М-т^, ZhUk"=T1Ih3', Zkk= —ютл^,
Z (10)_т I (5) Z (11)_I (6) т(7) Z (12) I (8)
^hk ~ 1 21 hk > ^hk ~xhk 1 2j^hk 1hk
ZT_- ЮТ2 Ihf, Zh221)_0, Zh222)_Ih6)-T2 If^,
Qh'^)_ 0, если l ^ j, 0 < l, j < 2;
q(00)_ г(2) q(11) _ Q(22)_ г(7)
hk hk hk hk hk
k_0
k_0
Здесь 0 < ) < 9, обозначают
= фФ1)+2а\ф'к, ф' )+а4(Фк, ф„), ¡Ц =- ¡а[(и 'ф', фк )+(и' ф, фк) +
+(и Фк, Ф )+а2(и фк, фк)], =-ф,ф')-а2ф,ф„), ¡(3)=-(^к,фк),
V), ¡¿в8)=-(иV'V), ¡^е
где т1=р/т, т2=1/т. Здесь
V' =- (2к + В)!'+1, фк =- (2к + 3) V' ф^=- (2к + 3)(2к + 5)!к+2.
(8)
Очевидно, что с помощью формул (8) и результатов работы [14] интегралы ¡кк могут быть рассчитаны аналитически. Например, первые два интеграла
в формуле для ¡'о) запишутся следующим образом: (фк", = 4(2к + 3)(2к + 3)(2к + 5)ёкк,
(фк, фк) = (2к + 3) 2(2к + 3)
2 +
2(2к + 3)
2к + 7
^кк -
2к + 7 6к'к-2 26к'к+2 В таблице для сравнения приведены значения собственного числа А0, отвечающего наиболее опасной моде, для уравнения Орра-Зоммерфельда в случае течения Пуазейля при а = 1, Яе = 104 (см. [15-17]).
Заключение. На рисунке 1 представлена часть точечного спектра для случая течения Куэтта-Пуазейля смеси при А = 0.1, а = 1, Яе = 104, Б = 10-5 и р = 0.1. Видно, что при малых величинах т наличие частиц в потоке проявляется лишь в малом смещении ограниченных частей спектра, которое может быть описано теорией аналитических возмущений линейных операторов. На рисунке 2 проиллюстрировано отличие спектра для монодисперсной смеси от случая однородной жидкости. Внутри области, отмеченной овалом, содержится бесконечное множество точек спектра задачи (3), о чем свидетельствует расходимость приближения.
На рисунке 3 показаны достаточные условия устойчивости для различных чисел Рейнольдса. Здесь, следуя традиции [12], через У обозначена величина
У = ^{А}/ а. Видно, что в отличие от однородной жидкости наблюдаются зоны метастабильности стационарного течения в закритической области [8, 9]. Необходимо отметить, что устойчивость основного течения смеси определяется дискретными собственными значениями.
0
К.е(Я)
-0.2
-0.4-
-0.6-
-0.8-
э _____9..... 9 О 3 С*. а У ф Ц]
С* и ^ о- \ укГ^
6 I 1т (Я)
0 0.2 0.4 0.6 0.8 1 Рис. 1. Спектр течения Куэтта — Пуазейля смеси при А = 0.1, а = 1, Яе = 104, 5 = 10-5 и р = 0.1
0 0.2 0.4 Гт(Я) 0.6 0.8 1
Рис. 2. Спектр течения Пуазейля смеси при а = 1, Яе = 104, 5 = 10-4 и р = 0.1
Ад = 0.23752649 +0.00373967/ Orszag [13]
А = 0.23752708 +0.00373980/ Бо^агга [13]
А = 0.2375264888204 +0.0037396706229/ Рор (К = 96, N = 512) [14]
А = 0.2375264888206 +0.0037396706230/ N = 64 [15]
Ад = 0.2375264888204 +0.00373967062291 N = 64
Рис. 3. Зависимость max;К(А/ a), Y = ÎK(A / а), от числа Рейнольдса. S = : 0-чистая жидкость;
а
1-1. E-6; 2-3. E-6; 3-5. E-6; 4-1.5E-5; 5-2. E-5; 6-5. E-5; 7-1. E-4; 8-1. E-3
Напряжения Рейнольдса в континууме дисперсных частиц не меняют знак в окрестности критического слоя. На границе канала не возникает разность фаз колебаний у продольной и поперечной компонент возмущения поля скоростей, поскольку в облаке дисперсных частиц отсутствует вязкое трение [12]. При этом наблюдается картина, когда при фиксированных параметрах задачи может наблюдаться несколько чисел Рейнольдса К.е4, соответствующих ней-
тральной устойчивости (см. рис. 3). Эффект «окон» устойчивости в закритической области, т. е. зон мета-стабильности основного течения, наблюдался в широких диапазонах значений параметров течения [8]. Отметим результаты работ [18, 19], в которых рассмотрена задача Рэлея-Бенара для дисперсной смеси и обнаружено, что в областях метастабильности наряду с устойчивой неподвижной точкой может наблюдаться неустойчивое стохастическое множество.
Библиографический список
1. Crowe C.C.T., Schwarzkopf J.D., Sommerfeld M. Multiphase Flows with Droplets and Particles. CRC Press, Boca Raton, 2011. DOI: 10.1201/b11103
2. Peter J. Schmid Nonmodal Stability Theory // Annual Review of Fluid Mechanics. — 2007. — Vol. 39. DOI: 10.1146/annurev.fluid.38.050304.092139
3. Divakar Viswanath. Spectral integration of linear boundary value problems // Journal of Computational and Applied Mathematics. — 2015. — Vol. 290. D0I:10.1016/j.cam. 2015.04.043
4. Payel Dasa, Gnaneshwar Nelakantia, Guangqing Long. Discrete Legendre spectral projection methods for Fredholm — Hammerstein integral equations // Journal of Computational and Applied Mathematics. — 2015. — Vol. 278. DOI: 10.1016/j.cam.2014.10.012
5. Balachandar S. and John K. Eaton Turbulent Dispersed Multiphase Flow // Annual Review of Fluid Mechanics. — 2010. — Vol. 42.
6. Saber A., Lundstrom T.S. and Hellstrom J.G.I. Turbulent Modulation in Particulate Flow: A Review of Critical Variables // Engineering. — 2015. — Vol. 7.
7. Боронин С.А. Устойчивость плоского течения Куэт-та дисперсной среды с конечной объемной долей частиц // Известия РАН. МЖГ. — 2011. — № 1.
8. Никитенко Н.Г., Сагалаков А.М., Попов Д.И. О достаточных условиях устойчивости течения Куэтта — Пуа-зейля монодисперсной смеси // Теплофизика и аэромеханика. — 2011. — № 2.
9. Попов Д.И., Утемесов Р.М., Оценка правой границы спектра в задаче об устойчивости параллельного течения двухфазной жидкости // Известия Алтайского гос. ун-та. — 2012. — № 1-2 (73).
10. Гершуни Г.З., Жуховицкий Е.М., Непомнящий А.А. Устойчивость конвективных течений. — М., 1989.
11. Foias C., Manley O., Rosa R., Temam R. Navier — Stokes Equations and Turbulence. — Cambridge, 2004.
12. Drazin P.G., Reid W.R. Hydrodynamic stability. — Cambridge, 1981.
13. Canuto C., Quarteroni A. Approximation results for orthogonal polynomials in Sobolev spaces, Math. Comp. — 1982. — Vol. 38.
14. Shen J. Efficient Spectral-Galerkin Method I. Direct Solvers for the Second and Fourth Order Equations Using Legendre Polynomials // SIAM J. Sci. Comput. — 1994. — Vol. 15, № 6.
15. Dongarra J.J., Straughan B., Walker D.W., Chebyshev tau-QZ algorithm methods for calculating spectra of hydrodynamic stability problems, Appl. — 1996. — Numer. Math. 22.
16. Pop I.S., Gheorghiu C.I., A Chebyshev-Galerkin method for fourth order problems, in Approximation and Optimization, Proceedings of the International Conference on Approxi-
mation and Optimization, D.D. Stancu et al. (eds.), Transilva-nia Press, Cluj-Napoca. — 1997. — Vol. II.
17. Попов Д.И., Утемесов Р.М. Спектральные коэффициенты операторов дифференцирования высокого порядка для ограниченных и неограниченных интервалов // Известия Алтайского гос. ун-та. — 2011. — № 1-2 (69).
18. Попов Д.И., Утемесов Р.М., Маломодовое приближение в задаче Бенара-Рэлея для двухфазной смеси // Известия Алтайского гос. ун-та. — 2014. — № 1-1 (81). DOI:10.14258/izvasu(2014)1.1-51
19. Попов Д.И., Утемесов Р.М. Моделирование конвекции Бенара-Рэлея дисперсной смеси // Ломоносовские чтения на Алтае: фундаментальные проблемы науки и образования : сб. науч. ст. Междунар. конф., Барнаул, 20-24 октября, 2015. — Барнаул, 2015.
