Известия РАН. Механика жидкости и газа, 2020, № 5, стр. 33-45

НЕСТАЦИОНАРНЫЕ ДВУХФАЗНЫЕ ТЕЧЕНИЯ ГАЗА С ЧАСТИЦАМИ В РЕШЕТКАХ ПРОФИЛЕЙ

Д. А. Романюк a*, Ю. М. Циркунов a**

a Балтийский государственный технический университет “ВОЕНМЕХ” им. Д.Ф. Устинова
Санкт-Петербург, Россия

* E-mail: romanyuk-da@rambler.ru
** E-mail: tsrknv@bstu.spb.su

Поступила в редакцию 16.01.2020
После доработки 08.02.2020
Принята к публикации 12.03.2020

Полный текст (PDF)

Аннотация

Численно исследуются нестационарные течения газа с твердыми частицами в системе двух плоских решеток “ротор–статор”. Течение несущего газа моделируется на основе уравнений Навье–Стокса (псевдо-DNS подход) и для сравнения на основе уравнений Рейнольдса (URANS подход) с k – ω SST моделью турбулентности Ментера. В обоих случаях уравнения решаются численно с помощью конечно-объемного метода второго порядка. Течение бесстолкновительной примеси моделируется методом Лагранжа, а столкновительной – методом Монте-Карло. Учитывается обратное влияние частиц на течение газовой фазы. Анализируется влияние различных факторов случайной природы (разброс частиц по размерам, рассеяние при отскоке от лопаток, столкновения между частицами) на картину течения примеси и профили ее концентрации на выходе из статорной решетки.

Ключевые слова: двухфазное течение газа с частицами, метод Лагранжа, метод Монте-Карло, разброс частиц по размерам, рассеяние частиц при отскоке, столкновения между частицами, численное моделирование

1. ВВЕДЕНИЕ

При полете летательного аппарата (ЛА) с газотурбинным двигателем в запыленной атмосфере дисперсные частицы могут засасываться в воздухозаборник и вызывать эрозию лопаток компрессора и последующих элементов тракта. Моделированию течения примеси, ее взаимодействию с лопатками компрессора и турбины, прогнозированию ресурса двигателя из-за эрозионного воздействия примеси в последнее время уделяется большое внимание [15]. Значительную опасность для ЛА представляет вынужденный полет в рассеянных облаках вулканического пепла [6]. Предсказание поведения частиц при обтекании ЛА и при течении в компрессоре и последующих элементах проточного тракта газотурбинного двигателя требует разработки адекватных математических моделей двухфазных течений. Наибольшие трудности здесь связаны с моделированием эффектов случайной природы, к которым относятся столкновения между частицами, рассеяние частиц при отскоке от обтекаемой поверхности из-за их несферической формы, разброс частиц по размерам, случайное положение частиц в невозмущенном потоке. Изучению этих эффектов посвящено очень небольшое число работ. В [7] была предложена кинетическая модель, описывающая динамику столкновительной примеси в известном поле течения несущего газа. В дальнейшем была разработана континуально-кинетическая модель газовзвеси [8], в которой несущий газ рассматривался как континуум, а примесь как дискретная среда. В этой модели учитывались столкновения между частицами и обратное воздействие примеси на течение несущего газа. Ударное взаимодействие частиц несферической формы со стенкой рассматривалось в [911]. Важным случайным фактором является разброс частиц по размерам [12, 13]. Все перечисленные факторы существенно влияют на картину течения примеси и, как следствие, на воздействие двухфазного потока на обтекаемую поверхность.

Течение в ступенях компрессора является, строго говоря, трехмерным, однако многие особенности течения могут быть изучены в рамках двумерной постановки. Такой подход оказался очень плодотворным и привел к созданию теории плоских решеток [14], которая успешно применяется для решения широкого круга задач в газодинамике турбомашин [15].

Целью данной работы является исследование роли случайных эффектов на распределение частиц в потоке в системе двух решеток “ротор–статор”.

2. ПОСТАНОВКА ЗАДАЧИ

Рассматривается плоское двухфазное течение газовзвеси через систему из двух решеток профилей (лопаток), которая в некотором смысле моделирует входную ступень осевого компрессора авиадвигателя. Первая решетка (ротор) движется с постоянной скоростью ${{V}_{r}}$ перпендикулярно невозмущенному потоку, а вторая решетка (статор) неподвижна (рис. 1). Обе решетки имеют одинаковый шаг s. Лопатки ротора установлены под углом $\beta $. Для анализа картины течения дисперсной фазы исследуется поведение облака частиц конечной длины h, в котором положение частиц задается случайным образом в соответствии с равномерным законом распределения.

Рис. 1.

Схема течения.

Среда частиц в облаке считается достаточно разреженной, так что можно пренебречь объемом дисперсной фазы в газовзвеси. Частицы являются более инерционными по сравнению с несущим газом, они не следуют линиям тока, могут сталкиваться с лопатками и отскакивать от них. Траектории отскочивших (отраженных) частиц могут пересекаться между собой и с траекториями падающих частиц. При увеличении концентрации примеси в невозмущенном потоке становятся существенными столкновения между падающими и отраженными частицами, что приводит к хаотизации движения примеси в пристеночном слое. При дальнейшем увеличении концентрации примеси она начинает влиять на течение несущего газа. Оценки концентрации частиц, когда нужно учитывать указанные эффекты, получены в [16]. При обтекании затупленных тел, каковыми являются лопатки в рассматриваемом течении, заметная роль столкновений между частицами проявляется при существенно более низкой концентрации примеси, чем та, при которой становится существенным обратное влияние примеси на течение несущего газа.

Динамическая чувствительность частиц к тонкой структуре течения несущего газа зависит от их относительной инерционности, которая характеризуется числом Стокса Stk = ${{\rho }_{p}}d_{p}^{2}{{V}_{\infty }}{\text{/}}(18{{\mu }_{\infty }}b)$. Этот параметр представляет собой отношение длины динамической релаксации частицы (при законе Стокса для силы сопротивления частицы) к характерному размеру течения. В данной работе в качестве последнего принята хорда лопаток $b$. Для частиц, рассмотренных в расчетах, число Стокса изменялось от нескольких единиц до нескольких десятков, что соответствует умеренно и высоко инерционным частицам. Это означает, что мелкомасштабная турбулентная структура течения несущего газа будет слабо влиять на поведение примеси.

3. МОДЕЛИРОВАНИЕ ТЕЧЕНИЯ НЕСУЩЕГО ГАЗА

В рассматриваемой задаче число Рейнольдса Re$_{\infty }$ = ${{\rho }_{\infty }}{{V}_{\infty }}b{\text{/}}{{\mu }_{\infty }} \sim {{10}^{6}}$ и, следовательно, течение является турбулентным. Традиционный подход в этом случае заключается в использовании осредненных уравнений Навье–Стокса (URANS- или комбинированный URANS/LES-подход) [17]. Оба подхода дают поля осредненных параметров газа. В то же время поведение примеси в потоке определяется, как известно, не осредненными, а актуальными значениями параметров несущей среды. Поэтому описание течения несущего газа на основе уравнений Навье–Стокса (DNS-подход) является в задачах о течении газовзвесей физически более корректным и, следовательно, предпочтительным для моделирования течения примеси. Строго говоря, численное моделирование турбулентных течений в этом случае требует нестационарной трехмерной постановки задачи и очень мелкой сетки (линейный размер ячеек должен быть порядка колмогоровского масштаба турбулентности) в расчетной области. В настоящее время такой подход реализуем только в простейших случаях. Для моделирования нестационарного течения несущего газа в системе решеток в данной работе использовался упрощенный подход: численно решались двумерные нестационарные уравнения Навье–Стокса (псевдо-DNS подход). Для сравнения было также выполнено решение нестационарных двумерных уравнений Рейнольдса (URANS подход) с k – ω SST моделью турбулентности Ментера [18, 19].

Уравнения Навье–Стокса для сжимаемого газа в рассматриваемом случае могут быть записаны в декартовых координатах $(x,y)$ в следующей векторной форме

$\frac{{\partial {\mathbf{Q}}}}{{\partial t}} + \frac{{\partial {{{\mathbf{F}}}_{x}}}}{{\partial x}} + \frac{{\partial {{{\mathbf{F}}}_{y}}}}{{\partial y}} = \frac{{\partial {{{\mathbf{G}}}_{x}}}}{{\partial x}} + \frac{{\partial {{{\mathbf{G}}}_{y}}}}{{\partial y}} + {\mathbf{H}}$
где Q – вектор консервативных переменных, Fx и ${{{\mathbf{F}}}_{y}}$ – векторы конвективных слагаемых, Gx и ${{{\mathbf{G}}}_{y}}$ – векторы диффузионных слагаемых. Они определяются соотношениями

${\mathbf{Q}} = \left( \begin{gathered} \rho \\ \rho {{{v}}_{x}} \\ \rho {{{v}}_{y}} \\ \rho e \\ \end{gathered} \right),\quad {{{\mathbf{F}}}_{x}} = \left( \begin{gathered} \rho {{{v}}_{x}} \\ \rho {v}_{x}^{2} + p \\ \rho {{{v}}_{x}}{{{v}}_{y}} \\ (\rho e + p){{{v}}_{x}} \\ \end{gathered} \right),\quad {{{\mathbf{F}}}_{y}} = \left( \begin{gathered} \rho {{{v}}_{y}} \\ \rho {{{v}}_{x}}{{{v}}_{y}} \\ \rho {v}_{y}^{2} + p \\ (\rho e + p){{{v}}_{y}} \\ \end{gathered} \right)$
${{{\mathbf{G}}}_{x}} = \left( \begin{gathered} 0 \\ {{\tau }_{{xx}}} \\ {{\tau }_{{xy}}} \\ {{{v}}_{x}}{{\tau }_{{xx}}} + {{{v}}_{y}}{{\tau }_{{xy}}} - {{q}_{x}} \\ \end{gathered} \right),\quad {{{\mathbf{G}}}_{y}} = \left( \begin{gathered} 0 \\ {{\tau }_{{xy}}} \\ {{\tau }_{{yy}}} \\ {{{v}}_{x}}{{\tau }_{{xy}}} + {{{v}}_{y}}{{\tau }_{{yy}}} - {{q}_{y}} \\ \end{gathered} \right)$

Здесь

${{\tau }_{{xx}}} = \frac{2}{3}\mu \left( {2\frac{{\partial {{{v}}_{x}}}}{{\partial x}} - \frac{{\partial {{{v}}_{y}}}}{{\partial y}}} \right),\quad {{\tau }_{{yy}}} = \frac{2}{3}\mu \left( {2\frac{{\partial {{{v}}_{y}}}}{{\partial y}} - \frac{{\partial {{{v}}_{x}}}}{{\partial x}}} \right)$
${{\tau }_{{xy}}} = {{\tau }_{{yx}}} = \mu \left( {\frac{{\partial {{{v}}_{x}}}}{{\partial y}} + \frac{{\partial {{{v}}_{y}}}}{{\partial x}}} \right)$
${{q}_{x}} = - \lambda \frac{{\partial T}}{{\partial x}},\quad {{q}_{y}} = - \lambda \frac{{\partial T}}{{\partial y}},\quad p = \rho RT,\quad e = {{c}_{V}}T + (1{\text{/}}2)({v}_{x}^{2} + {v}_{y}^{2})$

В записанных выше уравнениях $t$ – время; ${{{v}}_{x}}$ и ${{{v}}_{y}}$ – компоненты вектора скорости вдоль осей x и $y$; $\rho $, $p$, $e$, $T$, $\mu $ и $\lambda $ – плотность, давление, удельная полная энергия, температура, коэффициенты вязкости и теплопроводности газа соответственно; $R$ – газовая постоянная; ${{c}_{{v}}}$ и cp – удельные теплоемкости при постоянном объеме и давлении (для воздуха $R$ = 287 Дж/(кг · K), ${{c}_{p}}$ = = 1000 Дж/(кг · K); $\gamma = {{c}_{p}}{\text{/}}{{c}_{{v}}}$ = 1.4. Значение $\mu $ вычислялось по формуле Сазерленда: $\mu = {{\mu }_{s}}{{(T{\text{/}}{{T}_{s}})}^{{3/2}}}$(Ts + Cs)/(T + Cs), где ${{\mu }_{s}} = 1.71 \times {{10}^{{ - 5}}}$ Н · с/м2, ${{T}_{s}} = 273.15$ K, ${{C}_{s}} = 117.0$ K; $\lambda = {{c}_{p}}\mu {\text{/Pr}}$; число Прандтля было принято равным Pr = 0.71.

Моделирование источникого члена H, который описывает обратное влияние дисперсной примеси на несущий газ, и граничные условия будут описаны ниже.

Уравнения Рейнольдса и использованная в данной работе k – ω SST модель турбулентности подробно описаны в [18, 19], здесь они не приводятся ввиду их громоздкости.

4. МОДЕЛИРОВАНИЕ ТЕЧЕНИЯ ДИСПЕРСНОЙ ФАЗЫ

При моделировании движения примеси в потоке несущего газа все частицы рассматриваются как сферические. Такое допущение может быть использовано и в случае несферических частиц, если они выпуклые и изометрические, т.е. имеют близкие размеры по всем пространственным направлениям (см., например, [20]). В данной работе для оценки роли рассеяния частиц при отскоке от лопаток наряду со сферическими частицами рассматривались и несферические частицы в виде прямоугольных параллелепипедов с близкими сторонами по всем трем измерениям.

В случае бесстолкновительной примеси частицы движутся независимо друг от друга, и для описания их движения использовался классический лагранжев подход [21]. Уравнение импульса и момента импульса для отдельной частицы имеет вид

${{m}_{p}}\frac{{d{{{\mathbf{v}}}_{p}}}}{{dt}} = {{{\mathbf{f}}}_{D}} + {{{\mathbf{f}}}_{M}},\quad {{J}_{p}}\frac{{d{{{\mathbf{\omega }}}_{p}}}}{{dt}} = {{{\mathbf{l}}}_{p}},\quad \frac{{d{{{\mathbf{r}}}_{p}}}}{{dt}} = {{{\mathbf{v}}}_{p}}$
${{{\mathbf{f}}}_{D}} = \frac{1}{2}{{C}_{D}}\pi r_{p}^{2}\rho \left| {{\mathbf{v}} - {{{\mathbf{v}}}_{p}}} \right|({\mathbf{v}} - {{{\mathbf{v}}}_{p}})$
${{{\mathbf{f}}}_{M}} = \frac{4}{3}{{C}_{\omega }}\pi r_{p}^{3}\rho \left[ {({\mathbf{\omega }} - {{{\mathbf{\omega }}}_{p}}) \times ({\mathbf{v}} - {{{\mathbf{v}}}_{p}})} \right]$
${{{\mathbf{l}}}_{p}} = \frac{1}{2}{{C}_{L}}r_{p}^{5}\rho \left| {{\mathbf{\omega }} - {{{\mathbf{\omega }}}_{p}}} \right|({\mathbf{\omega }} - {{{\mathbf{\omega }}}_{p}}),\quad {\mathbf{\omega }} = (1{\text{/}}2){\text{rot}}\,{\mathbf{v}}$
где ${{m}_{p}} = (4{\text{/}}3){{\rho }_{p}}\pi r_{p}^{3}$, ${{J}_{p}} = (2{\text{/}}5){{m}_{p}}r_{p}^{2}$, ${{{\mathbf{v}}}_{p}}$, ${{{\mathbf{\omega }}}_{p}}$, ${{{\mathbf{r}}}_{p}}$ – масса, момент инерции, скорость, угловая скорость и радиус-вектор частицы, ρp – плотность вещества частицы.

В уравнении импульса учитываются сила аэродинамического сопротивления частицы fD и поперечная сила Магнуса fM, которая важна при поступательном движении вращающейся частицы относительно несущего газа. Сильное вращательное движение частиц возникает при касательных ударах о лопатки решеток и при столкновениях частиц друг с другом. В уравнении момента импульса учитывается аэродинамический момент ${{{\mathbf{l}}}_{p}}$.

Другие составляющие межфазной силы (сила Сэффмана, эффект присоединенной массы, наследственная сила Бассэ и др.) в рассматриваемых диапазонах параметров (числа Рейнольдса в поступательном и вращательном движении частиц, число Стокса, отношение плотностей газа и вещества частиц и др.), как показали предварительные оценки, несущественны и поэтому не учитываются.

Коэффициент CD вычисляется по аппроксимационной формуле Хендерсона [22] для режима дозвукового обтекания частицы. Эта формула справедлива в широких диапазонах чисел Маха и Рейнольдса частицы при ее движении относительно несущей среды. Они учитывают инерционность, сжимаемость и разреженность несущего газа и имеют вид

$\begin{gathered} {{C}_{D}}({{\operatorname{Re} }_{p}},{{{\text{M}}}_{p}},{{T}_{p}}{\text{/}}T) = \\ = \;24\mathop {\left\{ {{{{\operatorname{Re} }}_{p}} + \sqrt {\frac{\gamma }{2}} {{{\text{M}}}_{p}}\left[ {4.33 + \frac{{3.65 - 1.53{{T}_{p}}{\text{/}}T}}{{1 + 0.353{{T}_{p}}{\text{/}}T}}exp\left( { - 0.247\sqrt {\frac{2}{\gamma }} \frac{{{{{\operatorname{Re} }}_{p}}}}{{{{{\text{M}}}_{p}}}}} \right)} \right]} \right\}}\nolimits^{ - 1} + \\ \end{gathered} $
$ + \;\left[ {\frac{{4.5 + 0.38(0.03{\text{R}}{{{\text{e}}}_{p}} + 0.48\sqrt {{{{\operatorname{Re} }}_{p}}} )}}{{1 + 0.03{{{\operatorname{Re} }}_{p}} + 0.48\sqrt {{{{\operatorname{Re} }}_{p}}} }} + 0.1{\text{M}}_{p}^{2} + 0.2{\text{M}}_{p}^{8}} \right]exp\left( { - \frac{{{{{\text{M}}}_{p}}}}{{2\sqrt {{{{\operatorname{Re} }}_{p}}} }}} \right) + $
$ + \;0.6\sqrt {\frac{\gamma }{2}} {{{\text{M}}}_{p}}\left[ {1 - exp\left( { - \frac{{{{{\text{M}}}_{p}}}}{{{{{\operatorname{Re} }}_{p}}}}} \right)} \right]$
где $\mathop {{\text{Re}}}\nolimits_p = 2\rho \left| {{\mathbf{v}} - {{{\mathbf{v}}}_{p}}} \right|{{r}_{p}}{\text{/}}\mu $ и ${{{\text{M}}}_{p}} = {{\left| {{\mathbf{v}} - {{{\mathbf{v}}}_{p}}} \right|} \mathord{\left/ {\vphantom {{\left| {{\mathbf{v}} - {{{\mathbf{v}}}_{p}}} \right|} {\sqrt {\gamma RT} }}} \right. \kern-0em} {\sqrt {\gamma RT} }}$ – числа Рейнольдса и Маха при движении частицы относительно несущего газа. Отношение ${{T}_{p}}{\text{/}}T$ в данной задаче принималось равным единице.

Для вычисления силы Магнуса используется комбинация точного решения [23] и аппроксимационной зависимости [24]

${{C}_{\omega }} = \left\{ \begin{gathered} 3{\text{/}}4,\quad 2{{\gamma }_{\omega }} < 0.45 \hfill \\ (3{\text{/}}8){{\gamma }_{\omega }}[0.45 + (2{{\gamma }_{\omega }} - 0.45)exp( - 0.075\gamma _{\omega }^{{0.4}}\mathop {{\text{Re}}}\nolimits_p^{0.7} )],\quad 2{{\gamma }_{\omega }} \geqslant 0.45 \hfill \\ \end{gathered} \right.$
${{\gamma }_{\omega }} = {{\omega }_{p}}{{r}_{p}}{\text{/}}\left| {{\mathbf{v}} - {{{\mathbf{v}}}_{p}}} \right|$

Выражение для коэффициента CL взято из работы [25]

${{C}_{L}} = \frac{{{{C}_{{l1}}}}}{{\sqrt {{{{\operatorname{Re} }}_{{p\omega }}}} }} + \frac{{{{C}_{{l2}}}}}{{{{{\operatorname{Re} }}_{{p\omega }}}}}$
где ${{\operatorname{Re} }_{{p\omega }}} = \rho \left| {\omega - {{\omega }_{p}}} \right|r_{p}^{2}{\text{/}}\mu $ – число Рейнольдса вращательного движения частицы, а значения ${{C}_{{l1}}}$ и ${{C}_{{l2}}}$ приведены в табл. 1.

Таблица 1.

Значения коэффициентов ${{C}_{{l1}}}$ и ${{C}_{{l2}}}$

Repω 0–6 6–20 20–50 50–4 × 104
Cl1 0 5.32 6.44 6.45
Cl2 16π 37.2 32.2 32.1

Моделирование ударного взаимодействия частиц с лопатками решеток выполнено для сферических частиц и частиц в форме вытянутых прямоугольных призм с отношением каждой из двух коротких сторон к длинной, равным 0.8. При отскоке несферические частицы хаотически рассеиваются ввиду их случайной ориентации в момент удара.

Для описания ударного взаимодействия отдельной частицы с лопатками использовались уравнения импульса и момента импульса в интегральной форме. Эти уравнения связывают поступательную и угловую скорости частицы до и после удара. Для ее замыкания задавались коэффициенты восстановления нормальной ${{a}_{{pn}}} = - {v}_{{pn}}^{ + }{\text{/}}{v}_{{pn}}^{ - }$ и касательной ${{a}_{{p\tau }}} = {v}_{{p\tau }}^{ + }{\text{/}}{v}_{{p\tau }}^{ - }$ компонент вектора скорости (верхние индексы “–” и “+” относятся к параметрам частицы до и после столкновения). Для сферических частиц рассматривались коэффициенты восстановления скорости центра масс частицы, которые вычислялись по эмпирическим формулам [26, 27], учитывающим зависимость коэффициентов от скорости удара и угла падения (в качестве материала лопаток была взята сталь). Столкновения призматических частиц с лопатками рассматривались в трехмерной постановке и задавались коэффициенты восстановления компонент скорости не центра масс каждой частицы, а точки ее контакта при соударении [9]. Пространственная ориентация частицы перед столкновением принималась случайной и равновероятной. Удар считался нескользящим (коэффициент восстановления касательной скорости точки контакта задавался равным нулю). Коэффициент восстановления нормальной скорости точки контакта, описывающий степень неупругости удара, был принят равным 0.5. Учитывались возможные повторные соударения призматических частиц в процессе отскока.

При моделировании течения полидисперсной примеси был принят логарифмически-нормальный закон распределения массы частиц по размерам в невозмущенном потоке

${{g}_{\infty }}({{r}_{p}}) = \frac{1}{{\sqrt {2\pi } {{r}_{p}}ln\sigma }}exp\left[ { - {{{\left( {\frac{{ln{{r}_{p}} - ln{{r}_{g}}}}{{\sqrt 2 ln\sigma }}} \right)}}^{2}}} \right]$

Здесь ${{r}_{g}} = {{r}_{{pg}}}exp(l{{n}^{2}}\sigma )$, lnσ – среднеквадратичное отклонение логарифма радиуса частиц, ${{r}_{{pg}}}$ – наиболее вероятный радиус частиц.

Отраженные от поверхности лопаток частицы могут сталкиваться с падающими, что приводит к хаотизации движения всех вовлеченных частиц.

Для описания столкновительной примеси используется кинетическая модель, предложенная в работе [7]. Она основана на рассмотрении несущего газа как континуума, а дисперсной фазы как дискретного множества частиц. Поведение дисперсной фазы описывается обобщенным кинетическим уравнением Больцмана, учитывающим взаимодействие частиц с несущим газом, неупругие столкновения и трение между частицами.

5. ЧИСЛЕННЫЙ МЕТОД

Для численного решения уравнений несущего газа расчетная область разбивалась на два вычислительных блока (рис. 2). Левый блок двигался вместе с роторной лопаткой со скоростью Vr (рис. 1). Правый блок был неподвижен. На верхней и нижней границах каждого блока ставились условия периодичности решения для газа и для частиц по пространству, что возможно ввиду одинакового шага обеих решеток при отсутствии больших отрывных областей, распространяющихся на несколько шагов решетки. На границе между блоками использовалась специальная процедура согласования решений, основанная на переинтерполяции газодинамических параметров.

Рис. 2.

Вид расчетной области.

В каждом из блоков расчетной области непосредственно около профилей лопаток строилась сетка с четырехугольными ячейками, вытянутыми вдоль контура профиля и числом ячеек поперек пограничного слоя около 20. В остальной части расчетной области в каждом блоке генерировалась неструктурированная сетка с треугольными ячейками. Размеры ячеек обеих сеток на их общей границе в каждом из блоков были согласованы между собой. Общее количество ячеек в расчетной области в приведенных далее результатах расчетов составляло около 250 тысяч. Предварительно была исследована сходимость результатов при измельчении сетки и установлено, что увеличение количества ячеек практически не влияет на поля параметров. Численный алгоритм решения уравнений несущего газа был основан на явном методе конечного объема с использованием схемы Роу [28] для расчета невязких потоков на гранях ячеек и центральной конечно-разностной аппроксимации вязких членов. По временной координате использовался метод предиктор-корректор второго порядка аппроксимации.

На входной границе (перед роторной решеткой) задавались полная энтальпия h0 = γp/(γ – 1)${{\rho }_{\infty }} + V_{\infty }^{2}{\text{/}}2$ и энтропийная функция ${{p}_{\infty }}{\text{/}}\rho _{\infty }^{\gamma }$, компонента скорости ${{{v}}_{y}}$ принималась равной нулю. На выходной границе (за статорной решеткой) задавалось статическое давление, равное $1.2{{p}_{\infty }}$. На контурах профилей ставились условия равенства нулю нормальной и касательной компонент вектора скорости и постоянная температура ${{T}_{w}} = 288.15$ K. В качестве начальных условий задавались параметры невозмущенного потока во всей расчетной области. Характерным временем нестационарного течения является время, за которое роторная лопатка проходит один шаг решетки $\tau = s{\text{/}}{{V}_{r}}$. Выход течения газа на квазипериодический режим происходил при $\tau * = t{\text{/}}\tau $ ~ 20...25. После этого в расчетную область вводилось облако частиц, их скорость на входной границе принималась равной скорости газа.

Столкновения частиц между собой и с лопатками рассматривались в трехмерной постановке. Длина облака в невозмущенном потоке была равна $h = 5b$, ширина облака была принята равной шагу решеток s, поперечная толщина $\Delta z$ выбиралась из условия инвариантности решения относительно $\Delta z$ при дальнейшем увеличении толщины. Общее число моделирующих частиц в облаке зависело от их объемной концентрации ${{\alpha }_{\infty }}$ и $\Delta z$ и составляло 840 тыс. для частиц радиуса 5 и 10 мкм и 105 тыс. для частиц радиуса 20 мкм.

Источниковый член H в уравнениях Навье–Стокса, который моделирует обратное влияние дисперсной примеси на несущий газ, в каждой ячейке расчетной сетки определялся следующим образом [29]

${\mathbf{H}} = - \left\langle {{{{\mathbf{H}}}^{p}}} \right\rangle ,\quad {\mathbf{H}}_{i}^{p} = \left( {\begin{array}{*{20}{c}} 0 \\ {{{f}_{{px}}}} \\ {{{f}_{{py}}}} \\ {{{a}_{p}}} \end{array}} \right),\quad {\mathbf{H}}_{i}^{p} = - \sum\limits_{i = 1}^N {\frac{{{\mathbf{H}}_{i}^{p}}}{{{{\Omega }_{i}}}}} $
где $N$ – количество частиц в рассматриваемой ячейке расчетной сетки, ${\mathbf{H}}_{i}^{p}$ – вектор параметров для каждой отдельной i-й частицы, ${{f}_{{px}}}$ и ${{f}_{{py}}}$ – компоненты силы, действующей на частицу со стороны газа, ${{a}_{p}} = {{{\mathbf{f}}}_{p}} \cdot {{{\mathbf{v}}}_{p}} + {{{\mathbf{l}}}_{p}} \cdot {{\omega }_{p}} + {{q}_{p}}$ – скорость изменения полной энергии i-й частицы вследствие взаимодействия с газом, fp – сила аэродинамического сопротивления частицы, lp – момент, действующий со стороны газа на частицу, ${{{\mathbf{v}}}_{p}}$, ${{\omega }_{p}}$ – поступательная и угловая скорости частицы. В рассматриваемой задаче разница между температурами несущего газа и частиц очень мала, и тепловой поток ${{q}_{p}}$ от газа к частице принимался равным нулю. Суммирование параметров происходит по всем частицам в рассматриваемой ячейке в момент t.

Для численного моделирования течения примеси с учетом столкновений между частицами использовался алгоритм, основанный на методе Монте-Карло. Он подробно описан в [29]. Такой подход позволяет не проводить постоянную трассировку траекторий всех частиц, что существенно упрощает вычислительную процедуру и значительно сокращает время расчета. В случае, когда столкновения между частицами не учитываются, данный метод сводится к обычным траекторным расчетам частиц.

Для численного моделирования обеих фаз был разработан собственный программный код, реализованный в среде Intel Fortran.

6. РЕЗУЛЬТАТЫ РАСЧЕТОВ И ИХ АНАЛИЗ

Для лопаток и ротора и статора принят профиль NACA0012. Хорда профиля $b$ принята равной 100 мм, шаг обеих решеток s = 70 мм, угол установки β = 45°. Скорость набегающего потока ${{V}_{\infty }} = 200$ м/с, поперечная скорость решетки ротора ${{V}_{r}} = 150$ м/с. Несущий газ – воздух с плотностью ${{\rho }_{\infty }} = 1.23$ кг/м3, давление ${{p}_{\infty }} = {{10}^{5}}$ Па, материал частиц SiO2 с плотностью ${{\rho }_{p}} = 2650$ кг/м3. Параметры невозмущенного потока соответствуют числу Маха ${{{\text{M}}}_{\infty }} = 0.59$.

В расчетах радиус частиц монодисперсной примеси ${{r}_{p}}$ варьировался от 5 до 20 мкм, что соответствует числам Стокса Stk $ = 2{{\rho }_{p}}r_{p}^{2}{{V}_{\infty }}{\text{/}}(9{{\mu }_{\infty }}b)$ от 1.72 до 27.55. Для полидисперсных частиц параметр σ в логарифмически-нормальном законе распределения был принят равным 1.2, величина наиболее вероятного радиуса ${{r}_{{pg}}}$ варьировалась от 5 до 20 мкм.

Рассмотрим сначала результаты численного моделирования для несущего газа. На стадии тестирования программы было проведено сравнение полей параметров несущего газа, получающихся при численном решении уравнений Навье–Стокса и уравнений Эйлера (вязкость газа равна нулю). Оно показало, что во втором случае энтропийная функция $S = p{\text{/}}{{\rho }^{\gamma }}$, которая чувствительна к диссипативным эффектам, практически постоянна во всем поле течения, что говорит о малой схемной “вязкости” метода расчета. Последующее сравнение решений уравнений Навье–Стокса (псевдо-прямое численное моделирование) и уравнений Рейнольдса показало, что в первом случае получается значительно более подробная вихревая структура потока, в то время как при использовании уравнений Рейнольдса вихревая структура следов за лопатками полностью “размазана” (рис. 3). Поля энтропийной функции на рис. 3 соответствуют моменту t = 3 мс от начала счета, когда нестационарные процессы выходят на квазипериодический режим. Диапазон изменения энтропийной функции для уравнений Навье–Стокса получился около 10%, а для уравнений Рейнольдса 1.25%, что свидетельствует о значительно более высоком разрешении вихревой структуры потока в первом случае, так как именно в пограничных слоях и в вихревых следах за лопатками происходит диссипация, что и приводит к росту энтропии.

Рис. 3.

Поле энтропийной функции $S = p{\text{/}}{{\rho }^{\gamma }}$, ${{V}_{\infty }} = 200$ м/с, момент времени $t = 3$ мс: (а) – уравнения Навье–Стокса, (б) – уравнения Рейнольдса.

Наряду с полями энтропийной функции была построена зависимость числа Маха от времени в точке на расстоянии 2 мм перед профилем статорной решетки для обеих моделей течений несущего газа (рис. 4). В случае модели Навье–Стокса наблюдается хаотически-периодический характер изменения. Это связано с тем, что наряду с периодическим прохождением лопаток роторной решетки относительно лопаток статорной (период равен $\tau = s{\text{/}}{{V}_{r}}$) возникает еще один периодический процесс – срыв вихрей за лопатками и образование вихревой дорожки, которая по структуре подобна классической дорожке Кармана. Период этого процесса зависит от числа Рейнольдса при обтекании отдельной лопатки, и он не равен $\tau $. Интерференция двух этих периодических процессов приводит к сложному характеру изменения поля течения во времени и, как следствие, к появлению хаотической составляющей для зависимости числа Маха от времени. Для уравнений Рейнольдса эта зависимость имеет регулярный периодический характер (рис. 4б), что связано с полной диссипацией вихревой структуры следов под действием так называемой турбулентной вязкости, которая на несколько порядков превосходит физическую вязкость.

Рис. 4.

Зависимость от времени числа Маха в точке перед профилем статорной решетки на расстоянии 2 мм, скорость набегающего потока ${{V}_{\infty }} = 200$ м/с: (а) – уравнения Навье–Стокса, (б) – уравнения Рейнольдса.

Результаты численного моделирования течения дисперсной фазы, представленные далее, получены при использовании уравнений Навье–Стокса для несущего газа. Во всех случаях, если не оговорено другое, объемная концентрация примеси в невозмущенном потоке была равна ${{\alpha }_{\infty }} = {{10}^{{ - 4}}}$.

На рис. 5 приведены мгновенные картины распределения монодисперсных сферических частиц в системе решеток через 3 мс с момента начала движения облака через входную границу расчетной области после установления квазипериодического течения несущего газа. Там же показаны осредненные за интервал с 4 по 6 мс профили относительной объемной концентрации частиц в сечении AB (рис. 1), расположенном на расстоянии $l = 50$ мм за статорной решеткой.

Рис. 5.

Мгновенные картины распределения монодисперсных сферических частиц в потоке (слева) и осредненные по времени профили относительной объемной концентрации примеси в сечении AB (справа): (а) – ${{r}_{p}} = 5$ мкм (Stk = 1.72), (б) – ${{r}_{p}} = 10$ мкм (Stk = 6.89), (в) – ${{r}_{p}} = 20$ мкм (Stk = 27.55).

Структура течения примеси сильно зависит от размера частиц (от их инерционности, определяемой числом Стокса). Исходное однородное облако частиц при движении через систему решеток существенно деформируется: образуются узкие слои с высокой концентрацией примеси и области с пониженной концентрацией, что особенно выражено для частиц радиуса ${{r}_{p}} = 5$ и 10 мкм. Справа от картин распределения частиц в потоке показаны осредненные профили концентрации примеси, из которых видно, что в рассматриваемом поперечном сечении средняя концентрация частиц радиуса ${{r}_{p}} = 5$ и 10 мкм может превышать в 2 и более раз концентрацию примеси в невозмущенном потоке. Для частиц радиуса 20 мкм профиль более равномерный, что связано с более равномерным распределением в потоке отраженных от лопаток частиц. Отметим, что мгновенные профили концентрации существенно более неравномерные, чем осредненные: для примера на рис. 6 приведены мгновенные профили для момента t = 5 мс. Видно, что они имеют множество пиков, причем концентрация частиц всех размеров в отдельных точках превышает невозмущенное значение ${{\alpha }_{\infty }}$ в 4 раза.

Рис. 6.

Мгновенные профили относительной объемной концентрации частиц в сечении AB в момент $t = 5$ мс.

Мгновенные картины распределения частиц в системе решеток для того же момента времени с учетом различных случайных эффектов (разброс частиц по размерам, столкновения между частицами, рассеяние несферических частиц при отскоке от поверхности) приведены на рис. 7. Как видно из сравнения с рис. 5, все эти эффекты приводят к сглаживанию распределения примеси в системе решеток и, как следствие, профилей концентрации, однако их влияние различно. Столкновения между частицами радиуса 20 мкм качественно изменяют профиль, в то время как разброс частиц по размерам и их рассеяние при отскоке из-за несферической формы приводят лишь к некоторому количественному перераспределению концентрации. Для частиц меньших размеров (${{r}_{p}}$ и ${{r}_{{pg}}}$ 5 и 10 мкм) сглаживание пиков на осредненных профилях еще более выражено.

Рис. 7.

Мгновенные картины распределения частиц в потоке (слева) и осредненные по времени профили относительной объемной концентрации примеси на линии AB (справа): (а) – сферические полидисперсные частицы, наиболее вероятный радиус 20 мкм, ${{\alpha }_{\infty }} = {{10}^{{ - 4}}}$, параметр дисперсии σ = 1.2; (б) – сферические монодисперсные частицы радиуса 20 мкм, ${{\alpha }_{\infty }} = {{10}^{{ - 3}}}$, учитываются столкновения между частицами; (в) – призматические монодисперсные частицы, длина стороны 40 мкм, ${{\alpha }_{\infty }} = {{10}^{{ - 4}}}$.

Столкновения между частицами для всех рассмотренных размеров становятся заметными при их объемной концентрации в невозмущенном потоке ${{\alpha }_{\infty }} \sim {{10}^{{ - 3}}}$, последняя величина соответствует отношению массовых концентраций примеси и газовой фазы порядка единицы. Что касается учета обратного влияния примеси на течение несущего газа, то в исследованном диапазоне ${{\alpha }_{\infty }} = {{10}^{{ - 5}}}\, \ldots \;{{10}^{{ - 3}}}$ этот учет практически не сказывается на картинах распределения частиц.

Аналогичные расчеты движения примеси были выполнены для тех же исходных параметров течения, но для поля течения несущего газа, получаемого при решении уравнений Рейнольдса. Результаты показали, что для частиц радиуса 10 и 20 мкм картина течения примеси практически не изменялась, для частиц радиуса 5 мкм отличия были заметные, но очень слабые. Это связано с достаточно высокой инерционностью рассмотренных частиц, на движение которых сравнительно мелкомасштабная вихревая структура следов за лопатками не оказывала сколь-нибудь заметного влияния.

ЗАКЛЮЧЕНИЕ

Численное исследование нестационарного двумерного течения несущего газа в системе решеток лопаток “ротор–статор”, выполненное на основе полных уравнений Навье–Стокса (псевдо-прямое моделирование) и уравнений Рейнольдса с k – ω SST моделью турбулентности, показало, что в первом случае в возмущенной области получается значительно более тонкая вихревая структура потока, причем возмущенное течение даже при одинаковых шагах решеток не является строго периодическим. Причина этого заключается в интерференции двух периодических процессов – прохождение роторных лопаток относительно статорных (определяется только шагом решеток и скоростью роторной решетки) и образование вихревого следа за лопатками (определяется профилем лопаток и числом Рейнольдса). Во втором случае вихревые следы за лопатками полностью размазываются и течение получается строго периодическим.

Первоначально пространственно-однородное поле концентрации примеси в возмущенной области становится существенно неоднородным: в случае монодисперсных частиц образуются узкие слои с высокой концентрацией примеси (в 2–4 раза выше, чем в невозмущенном течении) и появляются области с концентрацией заметно ниже, чем в невозмущенном течении, или даже свободные от частиц. Учет эффектов случайной природы (разброс частиц по размерам, столкновения между частицами, рассеяние частиц при отскоке из-за их несферической формы) всегда приводит к размазыванию слоев с высокой концентрацией частиц и при некоторых параметрах примеси (размер частиц, концентрация в невозмущенном потоке) к их полному исчезновению и выравниванию профилей концентрации. Обратное влияние примеси на течение несущего газа в исследованных диапазонах размеров (${{r}_{p}}$ = 5–20 мкм) и концентрации частиц (${{\alpha }_{\infty }} = {{10}^{{ - 5}}}{\kern 1pt} - {\kern 1pt} {{10}^{{ - 3}}}$) оказалось несущественным.

Исследование выполнено при частичной финансовой поддержке РФФИ в рамках научного проекта № 20-08-00711.

Список литературы

  1. De Giorgi M.G., Campilongo S., Ficarella A. Predictions of operational degradation of the fan stage of an aircraft engine due to particulate ingestion // J. Eng. Gas Turbines Power. May 2015. V. 137. Issue 5. 052603, 15 p. https://doi.org/10.1115/1.4028769

  2. Cosher C.R., Dunn M.G. Comparison of the sensitivity to foreign particle ingestion of the GE-F101 and P/W-F100 engines to modern aircraft engines // J. Eng. Gas Turbines Power. December 2016, V. 138. Issue 12. 121201, 9 p. https://doi.org/10.1115/1.4034021

  3. Saxena S., Jothiprasad G., Bourassa C., Pritchard B. Numerical simulation of particulates in multistage axial compressors // J. Turbomach. March 2017. V. 139. Issue 3. 031013, 9 p. https://doi.org/10.1115/1.4034982

  4. Döring F., Staudacher S., Koch C., Weiβschuh M. Modeling particle deposition effects in aircraft engine compressors // J. Turbomach. May 2017. V. 139. Issue 5, 051003, 10 p. https://doi.org/10.1115/1.4035072

  5. Suman A., Morini M., Kurz R., Aldi N., Brun K., Pinelli M., Spina P.R. Estimation of the particle deposition on a subsonic axial compressor blade // J. Eng. Gas Turbines Power. 2017. V. 139. № 1, 012604, 10 p. https://doi.org/10.1115/1.4034209

  6. Dunn M.G. Operation of gas turbine engines in an environment contaminated with volcanic ash // J. Turbomach. September 2012. V. 134. Issue 5, 051001, 18 p. https://doi.org/10.1115/1.4006236

  7. Волков А.Н., Циркунов Ю.М. Кинетическая модель столкновительной примеси в запыленном газе и ее применение к расчету обтекания тел // Изв. РАН. МЖГ. 2000. № 3. С. 81–97.

  8. Volkov A.N., Tsirkunov Y.M., Oesterle B. Numerical simulation of a supersonic gas-solid flow over a blunt body: the role of inter-particle collisions and two-way coupling effects // Int. J. Multiphase Flow. 2005. V. 31. Issue 12. P. 1244–1275. https://doi.org/10.1016/j.ijmultiphaseflow.2005.07.002

  9. Панфилов С.В., Циркунов Ю.М. Рассеяние несферических частиц примеси при отскоке от гладкой и шероховатой поверхностей в высокоскоростном потоке газовзвеси // Прикладная механика и техническая физика. 2008. Т. 49. № 2. С. 79–88.

  10. Зотиков А.С., Лашков В.А. Коэффициент восстановления скорости при ударе абсолютно упругой частицы в форме эллипсоида вращения // Вестник СПбГУ. Сер. 1. 2014. Т. 1 (59). № 2. С. 245–253.

  11. Arboleda B.Q., Qadir Z., Sommerfeld M., Beatove S.L. Modelling the wall collisions of regular non-spherical particles and experimental validation // Proc. of the ASME, FEDSM’2014. 2014. Paper № 21610, 12 p. https://doi.org/10.1115/FEDSM2014-21610

  12. Моллесон Г.В., Стасенко А.Л. Особенности обтекания затупленного тела сверхзвуковой полидисперсной струей с закруткой отраженных частиц // Теплофизика высоких температур. 2011. Т. 49. № 1. С. 73–80.

  13. Ревизников Д.Л., Способин А.В., Сухарев Т.Ю. Численное моделирование обтекания затупленного тела сверхзвуковым полидисперсным потоком // Теплофизика высоких температур. 2017. Т. 55. № 3. С. 418–425.

  14. Степанов Г.Ю. Гидродинамика решеток турбомашин. М.: Физматлит, 1962.

  15. Gostelow J.P. Cascade aerodynamics // Oxford: Pergamon, 1984. (Рус. перев.: Гостелоу Дж. Аэродинамика решеток турбомашин. М.: Мир, 1987. 392 с.)

  16. Tsirkunov Yu.M. Gas-particle flows around bodies – key problems, modeling and numerical analysis // Proc. 4th Int. Conf. Multiphase Flow (Ed. E. Michaelides). 2001. Paper 607. 31 p.

  17. Любимов Д.А. Анализ турбулентных струйных и отрывных течений в элементах ТРД комбинированными RANS/LES-методами высокого разрешения // Дисс. д.ф.-м.н. Институт прикладной математики им. М.В. Келдыша РАН. 2014. 289 с.

  18. Menter F.R. Two-equation eddy-viscosity turbulence models for engineering applications // AIAA Journal. 1994. V. 32. № 8. P. 1598–1605. https://doi.org/10.2514/3.12149

  19. Menter F.R., Kuntz M., Langtry R. Ten years of industrial experience with the SST turbulence model // Proc. 4th Int. Symp. Turbulence, Heat and Mass Transfer (eds.: K. Hanjalic, Y. Nagano, and M. Tummers). Begell House Inc., West Redding, 2003. P. 625–632.

  20. Hölzer A., Sommerfeld M. New simple correlation formula for the drag coefficient of non-spherical particles // Powder Technol. 2008. V. 184. Issue 3. P. 361–365. https://doi.org/10.1016/j.powtec.2007.08.021

  21. Crowe C.T., Sommerfeld M., Tsuji Yu. Multiphase flow with droplets and particles // Boca Raton: CRC Press, 1998. 471 p.

  22. Henderson Ch.B. Drag coefficient of spheres in continuum and rarefied flows // AIAA Journal. 1976. V. 14. № 6. P. 707–708. https://doi.org/10.2514/3.61409

  23. Rubinow S.I., Keller J.B. The transverse force on a spinning sphere moving in a viscous fluid // J. Fluid Mech. 1961. V. 11. Issue 3. P. 447–459. https://doi.org/10.1017/S0022112061000640

  24. Oesterlé B., Bui Dinh T. Experiments on the lift of a spinning sphere in a range of intermediate Reynolds numbers // Experim. in Fluids. 1998. V. 25. P. 16–22. https://doi.org/10.1007/s003480050203

  25. Dennis S.R.C., Singh S.N., Ingham D.B. The steady flow due to a rotating sphere at low and moderate Reynolds numbers // J. Fluid Mech. 1980. V. 101. Issue 2. P. 257–279. https://doi.org/10.1017/S0022112080001656

  26. Лашков В.А. Об экспериментальном определении коэффициентов восстановления скорости частиц потока газовзвеси при ударе о поверхность // Инженерно-физический журнал. 1991. Т. 60. № 2. С. 197–203.

  27. Циркунов Ю.М., Панфилов С.В., Клычников М.Б. Полуэмпирическая модель ударного взаимодействия дисперсной частицы примеси с поверхностью, обтекаемой потоком газовзвеси // Инженерно-физический журнал. 1994. Т. 67. № 5–6. С. 397–386.

  28. Roe P.L. Approximate Riemann solvers, parameter vector and difference schemes // J. Comp. Phys. 1981. V. 101. Issue 2. P. 357–372. https://doi.org/10.1016/0021-9991(81)90128-5

  29. Volkov A.N., Tsirkunov Yu.M. CFD/Monte Carlo simulation of collision-dominated gas-particle flows over bodies // Proc. of ASME, FEDSM’2002. 2002. Paper 31222. 14 p. https://doi.org/10.1115/FEDSM2002-31222

Дополнительные материалы отсутствуют.