1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 59, № 10, 2019

Журнал вычислительной математики и математической физики, 2019, T. 59, № 10, стр. 1792-1802

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

А. Д. Савельев *

ВЦ ФИЦ ИУ РАН
119991 Москва, ул. Вавилова, 4, Россия

* E-mail: savel-cc09@yandex.ru

Поступила в редакцию 01.04.2019
После доработки 17.05.2019
Принята к публикации 10.06.2019

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

Аннотация

На основе компактных разностных схем 14-го порядка аппроксимации проводится численное моделирование дозвукового нестационарного обтекания вязким газом прямой и обратной ступенек на плоской поверхности при числе Маха набегающего потока 0.1. Исследуются характеристики зон отрыва пограничного слоя в диапазоне изменения числа Рейнольдса от 103 до 107. В расчетах используются нестационарные уравнения Навье–Стокса, дополненные при необходимости формулами дифференциальной однопараметрической SA-модели турбулентной вязкости. Определяются условия возникновения пульсирующих зон отрыва пограничного слоя и дается физическое обоснование последующему переходу к турбулентному режиму течения. Библ. 22. Фиг. 8.

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

1. ВВЕДЕНИЕ

Отрыв пограничного слоя на плохообтекаемых препятствиях является часто встречающимся в инженерной практике явлением. Несмотря на большой объем накопленных экспериментальных данных [1], [2], задача математического предсказания возникающих при этом силовых и акустических нагрузок не теряет своей актуальности. Исторически задача численного моделирования препятствия в виде расположенной по потоку ступеньки берет начало от задачи течения жидкости в канале с внезапным расширением [3]. За счет резкого перепада давления возникает отрыв пограничного слоя на верхней кромке с образованием зоны возвратно-циркуляционного течения, характеристики которой зависят от состояния пограничного слоя [4]–[8]. В случае ступеньки, расположенной против набегающего дозвукового потока, отрывная зона возникает не только перед ним, как при сверхзвуковом внешнем течении, но и на верхней поверхности за угловой кромкой [9]. Простота геометрической формы и компактность расчетной области, присущая задачам о течениях вязкой жидкости, делают данную конфигурацию поверхностей привлекательной для отладки различных моделей турбулентности [10].

В случае дозвуковых течений вязкого газа актуальность данной задачи определяется не аэродинамическими нагрузками на твердую поверхность, а возникающим при некоторых условиях акустическим излучением [8], [11], [12]. Источником данного излучения служат нестационарные процессы в сдвиговых слоях, усиливающиеся за счет наличия зон рециркуляции течения. Фрагментация зон отрыва пограничного слоя и снос оторвавшихся вихревых структур вниз по потоку ведут непосредственно к трансформации ламинарного течения в турбулентное. В этом плане определение условий, при которых течение в отрывных зонах становится нестационарным, близко проблеме неустойчивости ламинарного пограничного слоя и его восприимчивости к внешнему воздействию, изучению которой посвящено большое количество экспериментальных, теоретических и вычислительных работ (см., например, [13]–17]).

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

2. ИСХОДНЫЕ УРАВНЕНИЯ, ПОСТАНОВКА ЗАДАЧИ И МЕТОД РАСЧЕТА

Для решения задачи будем использовать двумерные нестационарные осредненные по Рейнольдсу уравнения Навье–Стокса. Обезразмеренные по параметрам набегающего потока и характерному линейному размеру, они имеют следующий вид:

$\frac{\partial }{{\partial t}}\rho + \frac{\partial }{{\partial {{x}_{j}}}}(\rho {{u}_{j}}) = 0,\quad \frac{\partial }{{\partial t}}(\rho {{u}_{i}}) + \frac{\partial }{{\partial {{x}_{j}}}}(\rho {{u}_{i}}{{u}_{j}} + {{\sigma }_{{ij}}}) = 0,$
$\frac{\partial }{{\partial t}}(\rho e) + \frac{\partial }{{\partial {{x}_{j}}}}(\rho e{{u}_{j}} + {{u}_{i}}{{\sigma }_{{ij}}} - {{q}_{j}}) = 0.$
Здесь $t$ и ${{x}_{j}}$ – время и декартовы координаты, $\rho $, ${{u}_{j}}$, $e = {{\gamma }^{{ - 1}}}h + 0.5u_{i}^{2}$ – плотность, компоненты вектора скорости и полная энергия, $h$ – энтальпия, $\gamma $ – отношение удельных теплоемкостей.

Составляющие тензора напряжений ${{\sigma }_{{ij}}}$ и тепловой поток представляются следующим образом:

${{\sigma }_{{ij}}} = {{\delta }_{{ij}}}\left( {p + \frac{2}{3}\rho k} \right) - 2(\mu + {{\mu }_{t}}){{\operatorname{Re} }^{{ - 1}}}\left( {{{S}_{{ij}}} - \frac{1}{3}\frac{\partial }{{\partial {{x}_{k}}}}({{u}_{k}}{{\delta }_{{ij}}})} \right),$
${{q}_{j}} = {{\operatorname{Re} }^{{ - 1}}}(\mu {{\Pr }^{{ - 1}}} + {{\mu }_{t}}\Pr _{t}^{{ - 1}})\frac{\partial }{{\partial {{x}_{j}}}}h,\quad {{S}_{{ij}}} = 0.5\left( {\frac{{\partial {{u}_{i}}}}{{\partial {{x}_{j}}}} + \frac{{\partial {{u}_{j}}}}{{\partial {{x}_{i}}}}} \right),$
где ${{S}_{{ij}}}$ – компоненты тензора скоростей деформации. Здесь $p$ – давление, $k$ – кинетическая энергия турбулентности, $\mu $ и ${{\mu }_{t}}$ – коэффициенты молекулярной и турбулентной вязкостей, $\Pr = 0.72$ и ${{\Pr }_{t}} = 0.9$ – молекулярное и турбулентное числа Прандтля, $\operatorname{Re} $ – число Рейнольдса. Система дополняется уравнением состояния $p = {{\gamma }^{{ - 1}}}\left( {\gamma - 1} \right)\rho h$ и зависимостью коэффициента молекулярной вязкости от энтальпии в виде формулы Сазерленда [18]. Для турбулентных режимов течения используется однопараметрическая дифференциальная модель [19], которая базируется на уравнении переноса для турбулентной вязкости $\tilde {\nu }$:

$\begin{gathered} \frac{\partial }{{\partial t}}\left( {\rho \tilde {\nu }} \right) + \frac{\partial }{{\partial {{x}_{i}}}}\left( {\rho {{u}_{i}}\tilde {\nu }} \right) = \frac{1}{{\sigma \operatorname{Re} }}\frac{\partial }{{\partial {{x}_{i}}}}\rho \left( {\nu + \tilde {\nu }} \right)\frac{{\partial{ \tilde {\nu }}}}{{\partial {{x}_{i}}}} + \\ + \;{{C}_{{lt}}}\left[ {\frac{{{{C}_{{b2}}}\rho }}{{\sigma \operatorname{Re} }}\frac{{\partial{ \tilde {\nu }}}}{{\partial {{x}_{i}}}}\frac{{\partial{ \tilde {\nu }}}}{{\partial {{x}_{i}}}} + {{C}_{{b1}}}(1 - {{f}_{{t2}}})\rho \tilde {S}\tilde {\nu } - \frac{1}{{\operatorname{Re} }}({{C}_{{w1}}}{{f}_{w}} - {{C}_{{b1}}}{{f}_{{t2}}}{{{\text{k}}}^{{ - 2}}})\rho {{{\left( {\frac{{\tilde {\nu }}}{d}} \right)}}^{2}}} \right]. \\ \end{gathered} $

Коэффициент турбулентной вязкости ${{\mu }_{t}}$ определяется по формуле

${{\mu }_{t}} = \rho \tilde {\nu }{{f}_{{v1}}}$,
где

${{f}_{{{v}1}}} = \frac{{{{\chi }^{3}}}}{{({{\chi }^{3}} + C_{{{v}1}}^{3})}},\quad \chi = \frac{{\tilde {\nu }}}{\nu },\quad \nu = \frac{\mu }{\rho }.$

Модель имеет следующие вспомогательные функции и константы:

${{C}_{{w1}}} = \frac{{{{C}_{{b1}}}}}{{{{{\text{k}}}^{2}}}} + \frac{{1 + {{C}_{{b2}}}}}{\sigma },\quad g = r + {{C}_{{w2}}}({{r}^{6}} - r),\quad {{f}_{w}} = g{{\left( {\frac{{1 + C_{{w3}}^{6}}}{{{{g}^{6}} + C_{{w3}}^{6}}}} \right)}^{{1/6}}},$
${{f}_{{t2}}} = {{C}_{{t3}}}\exp ( - {{C}_{{t4}}}{{\chi }^{2}}),\quad r = \frac{1}{{\operatorname{Re} }}\frac{{\rho \tilde {\nu }}}{{\tilde {S}{{{\left( {{\text{k}}d} \right)}}^{{\text{2}}}}}},$
$\tilde {S} = S + \frac{1}{{\operatorname{Re} }}\frac{{\rho \tilde {\nu }}}{{{{{\left( {{\text{k}}d} \right)}}^{2}}}}{{f}_{{{v}2}}},\quad S = \sqrt {2{{S}_{{ij}}}{{S}_{{ij}}}} ,\quad {{f}_{{{v}2}}} = 1 - \frac{\chi }{{1 + \chi {{f}_{{{v}1}}}}},$
${{C}_{{b1}}} = 0.1355,\quad {{C}_{{b2}}} = 0.622,\quad \sigma = 2{\text{/}}3,\quad {{C}_{{w2}}} = 0.3,\quad {{C}_{{w3}}} = 2,$
${{C}_{{{v}1}}} = 7.1,\quad {{C}_{{t3}}} = 1.2,\quad {{C}_{{t4}}} = 0.5,\quad k = 0.41,\quad 0 \leqslant {{C}_{{lt}}} \leqslant 1,$
где $d$ – расстояние по нормали от стенки.

Далее производится переход к обобщенной криволинейной системе координат для возможности расчета обтекания тел с искривленными и изломанными границами. Данная процедура подробно описана в [20]. Согласно ей область течения в физической плоскости $\left( {x,y} \right)$ отображается на единичный квадрат расчетной плоскости $\left( {\xi ,\eta } \right)$ с помощью преобразования общего вида $\xi = \xi \left( {x,y} \right)$, $\eta = \eta \left( {x,y} \right)$, $0 \leqslant \xi \leqslant 1$, $0 \leqslant \eta \leqslant 1$, а производные по декартовым координатам некой функции $f$ в уравнениях приобретают вид

$\frac{{\partial f}}{{\partial x}} = \frac{1}{{{{J}^{{ - 1}}}}}\left( {\frac{{\partial f{{y}_{\eta }}}}{{\partial \xi }} - \frac{{\partial f{{y}_{\xi }}}}{{\partial \eta }}} \right),\quad \frac{{\partial f}}{{\partial y}} = \frac{1}{{{{J}^{{ - 1}}}}}\left( {\frac{{\partial f{{x}_{\xi }}}}{{\partial \eta }} - \frac{{\partial f{{x}_{\eta }}}}{{\partial \xi }}} \right),$
где ${{J}^{{ - 1}}} = {{x}_{\xi }}{{y}_{\eta }} - {{x}_{\eta }}{{y}_{\xi }}$ – якобиан, а ${{x}_{\xi }}$, ${{y}_{\xi }}$, ${{x}_{\eta }}$ и ${{y}_{\eta }}$ – метрические коэффициенты преобразования. Таким образом, в уравнениях, использующихся для проведения расчетов, появляются метрические коэффициенты, а сами расчеты осуществляются на равномерной сетке, позволяющей применять компактные схемы высокого порядка.

Постановка задачи представлена на фиг. 1. Поток вязкого газа с числом Маха на бесконечности ${{{\text{M}}}_{\infty }} = 0.1$ обтекает конфигурацию твердых поверхностей, представляющую собой уступ, расположенный по потоку (схема $а$) или против потока (схема $б$). Начало системы координат расположено в передней критической точке $F$. За характерный размер $L$ принимается расстояние $FG$ от передней кромки пластины до ближайшей угловой точки. Высота уступа $h$ составляет $0.3L$.

Фиг. 1.

Постановка задачи.

Поскольку скорость набегающего потока существенно дозвуковая, внешние границы расчетной области отодвигаются на достаточно большие по сравнению с размерами препятствий расстояния. Так, левая граница $AB$ располагается на дистанции $30L$ от передней точки пластины $F$, что соответствует $100$ высотам уступа $GE$. Расстояния до верхней границы $BC$ и длина пластины за препятствием $ED$ равняются $50L$. На границе $AB$ задаются граничные условия на основе характеристик. На границах $BC$ и $CD$ ставятся условия вытекания. На твердой поверхности $FGED$-условия прилипания и нулевые значения коэффициента турбулентной вязкости, а на линии $AF$-симметрии потока в вертикальном направлении.

В расчетах ламинарного обтекания используется сетка с $301 \times 201$ сгущенными к твердой поверхности узлами. Минимальное расстояние между узлами составляет ${{10}^{{ - 4}}}L$. При расчетах турбулентных течений число узлов поперек пограничного слоя увеличивается до $251$, а минимальный шаг сетки уменьшается до $2.5 \times {{10}^{{ - 6}}}L$, что связано с необходимостью описания тонкого ламинарного подслоя.

Для расчета параметров потока в переходной области течения используется коэффициент ${{C}_{{lt}}}$, но при этом его значения не вычисляются, а задаются в зависимости от значения числа Re. В том случае, когда течение считается ламинарным, коэффициент ${{C}_{{lt}}}$, стоящий перед источниковыми членами уравнения турбулентной вязкости, равен 0. Уравнения модели турбулентности приобретают вид уравнений переноса и не оказывают сколь либо заметного влияния на решение. В случае турбулентного течения, на участке пластины $FE$ до значения координаты $x{\kern 1pt} ' = x{\text{/}}L = 0.3$, течение считается ламинарным, ${{C}_{{lt}}} = 0$. Переход задается при значениях параметра $\operatorname{Re} \,\left( {x{\kern 1pt} '\; - 0.3} \right)$ от $5 \times {{10}^{5}}$до ${{10}^{6}}$. При этом величина коэффициента ${{C}_{{lt}}}$ меняется линейно от 0 до 1. Ниже по потоку течение турбулентное, ${{C}_{{lt}}} = 1$. Значение турбулентной вязкости в набегающем потоке задается ${{\mu }_{{t\infty }}} = 0.1$, а уровень кинетической энергии турбулентности полагается равным $0$.

В расчетах используются схемы высокого порядка по пространственным переменным. В чем состоят отличия между разными схемами и чем обусловлена необходимость повышения их порядка аппроксимации? Дело в том, что конечно-разностные производные отличаются от их исходных дифференциальных аналогов. Так, разложение в ряд Тейлора разностного представления первой производной $\Delta _{1}^{{(n)}}f$ от некой гладкой функции $f$ по координате $x$ на сетке ${{\omega }_{h}} = \{ {{x}_{j}} = jh,\,h = {\text{const}}\} $ имеет вид

$\Delta _{1}^{{(n)}}f = \frac{{\partial f}}{{\partial x}} + \sum\limits_{k = n}^\infty {{{\alpha }_{k}}\frac{{{{\partial }^{{k + 1}}}f}}{{\partial {{x}^{{k + 1}}}}}} {{h}^{k}},$
где ${{\partial f} \mathord{\left/ {\vphantom {{\partial f} {\partial x}}} \right. \kern-0em} {\partial x}}$ – аппроксимируемая производная, $n$ – порядок аппроксимации разностной схемы, а коэффициенты $\left| {{{\alpha }_{k}}} \right| \to 0$ при $k \to \infty $. С увеличением порядка аппроксимации схемы $n$, члены под знаком суммы со значениями $k < n$ устраняются и получаемое не ее основе разностное решение приближается к решению исходных дифференциальных уравнений, поскольку влияние членов, остающихся под знаком суммы, снижается. Тут следует выделить два момента. Во-первых, остаточные члены с четными значениями степени $k$ при $h$ недиссипативны и оказывают влияние только на аппроксимацию производной. Чем выше значение $k$, тем точнее схема. Во-вторых, члены с нечетными значениями $k$ – это компоненты диффузного типа, ответственные за стабилизацию решения, элементы так называемой “схемной вязкости”. Их суммарное влияние в зависимости от значений коэффициентов ${{\alpha }_{k}}$ делает схему устойчивой или неустойчивой. Их задача – не дать рассчитываемым параметрам выйти за известные границы, о чем обычно говорят как о “развале” задачи. В то же время, излишне сильная диссипация вносит искажения в результаты численного моделирования, что не способствует изучению такого явления как отрыв пограничного слоя.

Для расчетов используются составные компактные схемы [21]. Данные схемы на сетке ${{\omega }_{h}} = \left\{ {{{x}_{j}} = jh,\,\,h = {\text{const}}} \right\}$ представляют первую производную функции $f$ по координате $x$ в виде суммы аппроксимирующей и стабилизирующей составляющих:

$f{\kern 1pt} ' = \left[ {\Delta _{1}^{{\left( l \right)}} + {{{\left( { - 1} \right)}}^{m}}s\frac{{\Delta _{2}^{m}}}{{{{c}_{{2m}}}}}} \right]\frac{f}{{2h}},$
где $l,m = 1,2\;...$ – коэффициенты для обозначения порядка центральной и повторной разностей, $s = \pm 1$ – параметр, учитывающий направление потока, а ${{C}_{{2m}}}$ – нормирующая константа, своя для каждого значения $m$. При этом операторами $\Delta _{1}^{{(l)}}$ и $\Delta _{2}^{m}$ обозначены первая производная $l$-го порядка и четная $2m$-я производная, описываемые с помощью оператора сдвига по сетке ${{T}_{n}}{{f}_{j}} = f({{x}_{j}} + nh)$:

$E = {{T}_{0}},\quad {{\Delta }_{1}} = {{T}_{1}} - {{T}_{{ - 1}}},\quad {{\Delta }_{2}} = {{T}_{1}} - 2E + {{T}_{{ - 1}}}.$

В данном случае в качестве аппроксимирующей составляющей используется центрально-разностный мультиоператор 14-го порядка [21]

$\Delta _{1}^{{\left( {14} \right)}} = \frac{1}{3}\sum\limits_{p = 1}^3 {{{{(1 + {{a}_{p}}{{\Delta }_{2}})}}^{{ - 1}}}} (1 + {{b}_{p}}{{\Delta }_{2}}){{\Delta }_{1}}$
при следующих значениях коэффициентов ${{a}_{p}}$ и ${{b}_{p}}$:

$\begin{gathered} {{a}_{1}} = {\text{0}}{\text{.2390484593463527,}}\quad {{b}_{1}} = - {\text{0}}{\text{.05674445887469773,}} \\ {{a}_{2}} = {\text{0}}{\text{.1624727238229966,}}\quad {{b}_{2}} = - {\text{0}}{\text{.03796675193891609,}} \\ {{a}_{3}} = {\text{0}}{\text{.06001727836911225,}}\quad {{b}_{3}} = {\text{0}}{\text{.02207959560705087}}{\text{.}} \\ \end{gathered} $

В качестве стабилизирующих операторов $\Delta _{2}^{m}$ используются 16-е производные расщепленных и ориентированных в соответствии с направлением характеристик векторов потоков. Расчет смешанных производных, производных в источниковых членах и метрических коэффициентов также проводится на основе симметричной схемы $\Delta _{1}^{{(14)}}$.

Помимо конвективных членов, в уравнениях динамики вязкого газа присутствуют диссипативные члены вида $\partial (\mu \partial f{\text{/}}\partial x){\text{/}}\partial x$. Обычно они описываются стандартным образом операторами второго порядка. Для турбулентных течений, при моделировании которых коэффициент вихревой вязкости меняется в широких пределах, становится актуальным более точное описание диффузных членов уравнений. В данном случае разностное представление дифференциального оператора на основе схем 14 -го порядка выглядит следующим образом:

${{\delta }^{{(14)}}}{{I}^{{(14)}}}\mu {{\delta }^{{(14)}}}f,$
где
${{\delta }^{{\left( {14} \right)}}} = \frac{1}{3}\sum\limits_{p = 1}^3 {{{{(1 + a_{p}^{\delta }{{\Delta }_{2}})}}^{{ - 1}}}} [(1 - 3b_{p}^{\delta })({{T}_{{1/2}}} - {{T}_{{ - 1/2}}}) + b_{p}^{\delta }({{T}_{{3/2}}} - {{T}_{{ - 3/2}}})],$
и
${{I}^{{\left( {14} \right)}}} = \frac{1}{3}3\sum\limits_{p = 1}^n {{{{(1 + a_{p}^{I}{{\Delta }_{2}})}}^{{ - 1}}}} [(1 - b_{p}^{I})({{T}_{{1/2}}} + {{T}_{{ - 1/2}}}) + b_{p}^{I}({{T}_{{3/2}}} + {{T}_{{ - 3/2}}})],$
формулы компактных разности и интерполяции в точках между узлами разностной сетки, коэффициенты которых имеют следующие значения:

$\begin{gathered} a_{1}^{\delta } = {\text{0}}{\text{.210947273633837529,}}\quad b_{1}^{\delta } = {\text{0}}{\text{.1839392747705121678,}} \\ a_{2}^{\delta } = {\text{0}}{\text{.118639542332554637,}}\quad b_{2}^{\delta } = {\text{0}}{\text{.05792412786637447797,}} \\ a_{3}^{\delta } = {\text{0}}{\text{.03117204101727738,}}\quad b_{3}^{\delta } = - {\text{0}}{\text{.0061045456532170918,}} \\ \end{gathered} $
$\begin{gathered} a_{1}^{I} = {\text{0}}{\text{.2376211084878023794,}}\quad b_{1}^{I} = {\text{0}}{\text{.0339458726411146286,}} \\ a_{2}^{I} = {\text{0}}{\text{.1528151167445392922,}}\quad b_{2}^{I} = {\text{0}}{\text{.02183073096350561515,}} \\ a_{3}^{I} = {\text{0}}{\text{.0470637747676583076,}}\quad b_{3}^{I} = {\text{0}}{\text{.00672339639537973544}}{\text{.}} \\ \end{gathered} $

Для интегрирования по времени используется неявный метод Гаусса $2$-го порядка по времени релаксации в линиях. Отношение временного и пространственного шагов в основной массе расчетов равняется $1$.

3. РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

Расчеты показали высокую чувствительность отрывного течения к уровню вязких сил в потоке (числу Рейнольдса).

На фиг. 2 и 3 представлены поля линии тока в виде линий равных значений функции тока

$\psi = \int {\rho u\partial y - \int {\rho v\partial x} } ,$
полученных для ступенек, расположенных по потоку и против потока соответственно. При построении изолиний использованы безразмерные координаты $x{\kern 1pt} ' = x{\text{/}}L$ и $y{\kern 1pt} ' = y{\text{/}}L$. Области замкнутых изолиний соответствуют зонам рециркуляции течения. Решения для значений числа $\operatorname{Re} = {{10}^{3}}$ отмечены литерой $а$, ${{10}^{4}}$ – $б$ и ${{10}^{6}}$ – $в$. Согласно полученным результатам, единые при низких ламинарных значениях $\operatorname{Re} $ зоны отрыва с увеличением числа Рейнольдса постепенно вытягиваются и теряют устойчивость. Далее происходит фрагментация зон рециркуляции с образованием вихревых структур, которые, постепенно диссипируя, сносятся спутным потоком к выходной границе. Использование модели турбулентной вязкости с одновременным увеличением числа Рейнольдса приводит к тому, что зоны отрыва пограничного слоя в каждом случае вновь собираются воедино. Данная метаморфоза объясняется влиянием вязких сил, вновь усилившихся после введения добавочной турбулентной вязкости.

Фиг. 2.

Обтекание ступеньки, расположенной по потоку. Изолинии функции тока.

Фиг. 3.

Обтекание ступеньки, расположенной против потока. Изолинии функции тока.

Возникающие пульсации потока отражаются и на распределении давления на твердой поверхности. На фиг. 4 представлены изменения по времени $t{\kern 1pt} ' = tu{}_{\infty }{{L}^{{ - 1}}}$ значений коэффициента давления ${{C}_{p}} = 2(p - {{p}_{\infty }}){\text{/}}{{\rho }_{\infty }}u_{\infty }^{2}$ в угловой точке уступа, расположенного против потока, при следующих значениях числа Рейнольдса: 1 – Re = 103, 2 – 104, 3 – 105, 4 – 106. Изменение коэффициента давления за уступом, расположенным по потоку, в точке с координатами $x{\kern 1pt} ' = 3,\;y{\kern 1pt} ' = - 0.3$ при тех же значениях Re приведены на фиг. 5. При ламинарном режиме течения рост Re приводит к усилению пульсаций потока. Амплитуда колебаний, возникающих при Re > 103, с ростом дефицита вязких сил только увеличивается. Подключение модели турбулентной вязкости снижает вариативность пульсаций давления и делает их более регулярными. При Re = 107 колебания практически устраняются, а размеры зон рециркуляции потока снижаются за счет роста вязких сил. Последний результат противоречит экспериментальным данным, приведенным в [1], где отмечена тенденция к увеличению размеров отрыва турбулентного пограничного слоя с ростом $\operatorname{Re} $, что говорит о необходимости совершенствования существующих моделей турбулентной вязкости.

Фиг. 4.

Изменение по времени коэффициента давления в верхней угловой точке уступа, расположенного против потока.

Фиг. 5.

Изменение по времени коэффициента давления на поверхности за уступом, расположенным по потоку.

Колебания потока и движение вихревых структур вдоль поверхности вызывают излучение звуковых волн. На фиг. 6 представлено мгновенное поле нестационарного избыточного давления $\Delta \bar {p} = (p - \bar {p}){\text{/}}\bar {p}$ ($\bar {p}$ – среднее по времени давление), полученное для обратного уступа при значении $\operatorname{Re} = 3 \times {{10}^{5}}$. Зоны низкого давления соответствуют вихревым структурам, движущимся вдоль поверхности. Их постепенная диссипация в поле течения объясняется, помимо проходящих физических процессов, еще и увеличением размеров сеточных узлов, что в свою очередь вызвано необходимостью использования значительных размеров расчетного поля при низких значениях числа Маха набегающего потока.

Фиг. 6.

Поле нестационарного избыточного давления.

Большой интерес представляют условия, при которых начинается фрагментация сформировавшейся у верхней угловой кромки зоны отрыва пограничного слоя. Расчеты показали, что для уступа, расположенного по потоку, это происходит при значении числа Рейнольдса больше $4.5 \times {{10}^{3}}$, а в случае уступа, расположенного против потока – при Re > $1.1 \times {{10}^{3}}$. Градиент давления во втором случае больше, чем в первом. Непосредственное использование критерия отрыва ламинарного пограничного слоя [22], который в безразмерном виде можно представить как

${{C}_{B}} = \frac{{\operatorname{Re} {{\delta }^{{*2}}}}}{2}\frac{{d{{C}_{p}}}}{{dx}}$
($\delta {\kern 1pt} {\text{*}}$ – толщина вытеснения пограничного слоя перед точкой отрыва), дает в случае уступа, расположенного по потоку, значение $ \approx $1.4, а в другом случае довольно близкие значения 1.6–1.8. Можно предположить, что на основе данного критерия возможно, помимо возникновения собственно отрыва пограничного слоя, также предсказывать начало фрагментации и турбулизации возвратно-циркуляционного течения, хотя бы для случаев с фиксированной точкой отрыва.

При высоких ламинарных значениях Re за уступами наблюдается интенсивное перемещение вихревых структур вниз по потоку. Но интерес представляют не только мгновенные значения параметров течения, но и средние на неком временном промежутке. На фиг. 7 представлены профили осредненной по времени горизонтальной компоненты вектора скорости $\bar {u}$ у выходной границы расчетного поля. Цифры 1 и 3 относятся к уступам, расположенным, соответственно, по потоку и против потока при числе $\operatorname{Re} = 3 \times {{10}^{5}}$. Такие же профили, полученные для случая турбулентного течения при $\operatorname{Re} = {{10}^{6}}$, отмечены цифрами 2 и 4. Полученные профили осредненного течения больше похожи на турбулентные, чем на ламинарные. Уровни поверхностного трения также ближе к турбулентным. А вот толщина пограничного слоя значительно отличается от значений, полученных с использованием модели турбулентной вязкости как для уступа, расположенного по потоку, так и против потока. Она или меньше в первом случае, или существенно выше во втором. В то же время для турбулентного режима течения толщина пограничного слоя в обоих случаях имеет близкие значения.

Фиг. 7.

Профили осредненной по времени горизонтальной компоненты вектора скорости.

При проведении расчетов в ряде случаев были зафиксированы мелкомасштабные высокочастотные колебания давления. На фиг. 8 представлено изменение по времени коэффициента давления в передней кромке пластины (точка $F$ на фиг. 1), полученное в расчете обтекания уступа, расположенного по потоку при числе Рейнольдса 103. Временной участок, отмеченный буквой $а$, получен при соотношении временного и пространственного шагов 1 и выводе значений давления на каждой 5-й итерации по времени, буквой $б$ – при выводе на каждой итерации, $в$ – при соотношении временного и пространственного шагов 0.3. Более подробно описать данные колебания на используемой сетке не представлялось возможным. Если перейти к размерным величинам и положить, например, что длина пластины перед уступом составляет 0.1 м, то в последнем случае временной шаг составит 10–7 сек. Колебания давления лежат в диапазоне $3{\kern 1pt} - {\kern 1pt} 5 \times {{10}^{{ - 6}}}$ от давления в набегающем потоке.

Фиг. 8.

Колебания коэффициента давления на передней кромке пластины.

Полученные результаты требуют объяснения и обоснования. Естественно считать, что процессы вихреобразования, возникающие при решении нестационарных уравнений Навье–Стокса при высоких ламинарных значениях $\operatorname{Re} $, являются проявлением дефицита вязких сил в поле течения. Дальнейшее увеличение числа Рейнольдса ведет к усилению градиентов давления, что в свою очередь должно приводить к развалу решения. В природе этого не происходит, поскольку здесь следует подключение иного, отличного от ламинарного механизма межмолекулярного взаимодействия. То, что для успешных расчетов распределений давления, трения и теплопередачи используются именно модели на основе концепции добавочной турбулентной вязкости, служит этому дополнительным подтверждением. В свою очередь, современные модели турбулентности требуют дальнейшего совершенствования для описания нестационарных вихревых течений во всех необходимых практике диапазонах изменения чисел Маха и Рейнольдса.

Внешнее сходство профилей осредненной по времени горизонтальной составляющей вектора скорости с турбулентными не является подтверждением перехода к турбулентному режиму течения. Более надежным следует считать выполнение критерия отрыва турбулентного пограничного слоя [22].

Представленные на фиг. 8 в зоне $а$ временные зависимости коэффициента давления на передней кромке пластины на первый взгляд выглядят как турбулентные пульсации потока. В то же время при более высоких значениях $\operatorname{Re} $, когда пограничный слой тоньше, амплитуда пульсаций давления снижается. Поскольку структура изменения давления по времени не проясняется при более подробном выводе информации или уменьшении шага по времени, следует понимать, что это вычислительный эффект, связанный с конечностью шагов по времени и разностной сетке, а также с имеющим место разрывом параметров потока на передней кромке пластины.

4. ЗАКЛЮЧЕНИЕ

Проведенные на основе схем 14 -го порядка расчеты обтекания прямой и обратной ступенек вязким дозвуковым потоком газа в широком диапазоне изменения числа Рейнольдса показали сильную зависимость характеристик отрывного течения от уровня вязких сил в набегающем потоке. Определены значения $\operatorname{Re} $, при которых возникает распад единых отрывных зон и начинается формирование сбегающих вниз по потоку вихревых структур. Показано, что данные нестационарные эффекты являются следствием дефицита вязких сил в потоке. Снижение уровня пульсаций потока и формирование вновь единых зон отрыва пограничного слоя при использовании модели турбулентной вязкости следует рассматривать как эффект усиления вязких сил в поле течения.

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

  1. Чжен П. Отрывные течения. Т. 1. М.: Мир, 1972. 300 с.

  2. Попович К.Ф., Лаврухин Г.Н. Аэродинамика реактивных сопел. Т. 2. Обтекание донных уступов потоком газа. М.: Физматлит. 2009. 303 с.

  3. Armaly B.F., Durst F., Pereira J.C.F., Schonung B. Experimental and theoretical investigation of backward-facing step flow // J. of Fluid Mechanics. 1983. V. 127. P. 473–496.

  4. Елизарова Т.Г., Калачинская И.С., Шеретов Ю.В., Шильников Е.В. Численное моделирование отрывных течений за обратным уступом // Прикл. матем. и информатика. 2003. № 14. С. 85–118.

  5. Elizarova T.G., Nikolskii P.N., Lengrand J.C. A new variant of subgrid dissipation for LES method and simulation of laminar-turbulent transition in subsonic gas flows // Proceedings of the second Symposium on Hybrid RANS-LES Methods. 2007. Corfu. Greece. 17–18 June 2007. P. 73–74.

  6. Елизарова Т.Г., Никольский П.Н. Численное моделирование ламинарно-турбулентного перехода в течении за обратным уступом // Вестн. Московского ун-та. Серия 3. Физика. Астрономия. 2007. № 4. С. 14–17.

  7. Батенко С.Р., Терехов В.И. Влияние динамической предыстории потока на аэродинамику ламинарного отрывного течения в канале за обратным прямоугольным уступом // Прикл. механ. и техн. физ. 2002. Т. 43. № 6. С. 84–92.

  8. Lasher W.C., Taulbee D.B. On the computation of turbulent backstep flow // Int. j. Heat and Fluid Flow. 1992. V. 13. № 1. P. 30–40.

  9. Addad Y., Laurence D., Talotte C., Jacob M.C. Large eddy simulation of a forward-backward facing step for acoustic source identification // Int. j. Heat and Fluid Flow. 2003. V. 24. P. 562–571.

  10. Shur M., Strelets M., Zaikov L., Gulyaev A., Kozlov V., Sekundov A. Comparative numerical testing of one- and two-equation turbulence models for flows with separation and reattachment // AIAA Paper No 95-0863. 1995. 31 p.

  11. Gloerfelt X., Bailly C., Juve D. Computation of the noise radiated by a subsonic cavity using direct simulation and acoustic analogy // AIAA Paper No 2001–2226. 2001. 12 p.

  12. Голубев А.Ю., Ефимцов Б.М. Особенности структуры полей пульсаций давления в окрестности выступов // Известия РАН. Механ. жидкости и газа. 2015. № 1. С. 55–66.

  13. Wundrow D.W., Goldstein M.E. Effect on a laminar boundary layer of small-amplitude streamwise vorticity in the upstream flow // J. Fluid Mech. 2001. V. 426. P. 229–262.

  14. Chiu W.K., Norton M.P. The receptivity of laminar boundary layer to leading edge vibration // J. of Sound and Vibration. 1990. V. 141. P. 143–173.

  15. Luchini P. Reynolds-number independent instability of the boundary layer over a flat surface: optimal perturbations // J. Fluid Mech. 2002. V. 404. P. 289–309.

  16. Bandadi Y., Sbaibi A. Subsonic boundary layer receptivity to free stream acoustic perturbations // 13ėme Congrės de Mėcanique. 11–14 avril 2017. Meknės. Maroc. 3 p.

  17. Heinrich R.A., Choudhari M., Kerschen E.J. A comparison of boundary layer receptivity mechanisms // AIAA Paper No 88-3758. 1988. 8 p.

  18. Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1987. 840 с.

  19. Spalart P.R., Allmaras S.R. A one-equation turbulence model for aerodynamic flows // AIAA Paper No 92-0439. 1992. 21 p.

  20. Steger J.L. Implicit finite-difference simulation of flow about arbitrary two-dimensional geometries // AIAA J. 1978. V. 16. P. 678–685.

  21. Савельев А.Д. О мультиоператорном представлении составных компактных схем // Ж. вычисл. матем. и матем. физ. 2014. Т. 54. № 10. С. 1580–1593.

  22. Бам-Зеликович Г.М. Расчет отрыва пограничного слоя // Изв. АН СССР. ОТН. 1954. № 12. С. 68–85.

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

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал