Известия РАН. Механика жидкости и газа, 2023, № 3, стр. 36-46

СЖИМАЕМЫЕ ТЕЧЕНИЯ С ОСЦИЛЛИРУЮЩИМИ ГРАНИЦАМИ

О. А. Логвинов *

МГУ им. М.В. Ломоносова, Механико-математический факультет
Москва, Россия

* E-mail: oleglogvinov@mail.ru

Поступила в редакцию 24.10.2022
После доработки 16.12.2022
Принята к публикации 16.12.2022

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

Аннотация

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

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

Во многих технологических природных и физиологических процессах распространены течения в каналах с подвижными границами. Различные течения с осциллирующими твердыми частями возникают в авиационной промышленности. Отдельное, давно и успешно развивающееся в нашей стране направление связано с вибрационным воздействием на деформируемый пласт, насыщенный фильтрующейся жидкостью. В работах [1, 2] получено решение для стационарных колебаний двухфазной системы “пласт–жидкость” и показана возможность резкого увеличения амплитуды колебаний (резонанса) на определенных частотах для конкретных граничных условий. Установлено, что для всех параметров резонанс соответствует колебаниям пористого скелета и жидкости в одной фазе. Показано, что резонансные частоты дают возможность создавать значительные амплитуды на большом расстоянии от источника колебаний.

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

В работе [3] рассматривалось течение несжимаемой вязкой жидкости в плоском канале или цилиндрической трубе с жесткими стенками, движущимися по своей нормали, так что ширина канала или радиус трубы меняются со временем по заданному периодическому закону. В случае малой амплитуды колебаний стенок ${{\varepsilon }}$ (отнесенной к ширине или радиусу канала) решение получалось при произвольных вибрационных числах Рейнольдса α (произведение числа Рейнольдса на число Струхаля) разложением в ряды по $\varepsilon $. Интересно, что при малых ${{\alpha }}$ средние скорости индивидуальных частиц направлены к стенкам, тогда как при больших ${{\alpha }}$ дрейф происходит, в основном, от стенок. Исключение составляют частицы, образующие тонкий пограничный слой.

Также в [3] изучалось движение боковых стенок канала, возникающее в результате их упругого отклика на периодические давления, приложенные как на концах канала, так и на самих стенках. Оказалось, что система жидкость – стенки обладает собственной частотой. А в случае, если эта частота оказывается кратной частоте осцилляций внешнего давления, возможно возникновение резонанса. Предложенная постановка имеет приложение к течению крови в коронарных артериях больших млекопитающих.

В работе [4] рассматривалась линейная устойчивость решения, полученного в [3], для случая плоского канала. Было показало, что увеличение вибрационного числа Рейнольдса для достаточно больших амплитуд колебаний стенок ведет к детерминированному хаосу по сценарию Фейгенбаума через каскад удвоения периода.

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

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

В работе [7] изучался массоперенос в плоском канале с вибрирующими упругими стенками. Решение строилось разложением по малому параметру ${{\varepsilon }}$. Рассматривались не только колебания в виде бегущих волн, но и более общие виды движений стенок. Выяснилось, что при больших вибрационных числах Рейнольдса стоячие волны дают максимальный расход, на порядок превосходящий расход в случае бегущих по упругим границам волн. Таким образом, было опровергнуто устоявшееся мнение, что бегущие волны – оптимальный механизм прокачки жидкости в канале с колеблющимися деформируемыми стенками.

В [8] экспериментально моделировались течения воздуха в осциллирующей цилиндрической трубе. Частоты осцилляций давления в “покоящейся” трубе совпадают с ее собственными частотами. При вынужденных же осцилляциях стенок соответствующие частоты осцилляций давления являются высшими гармониками для частот стенок. Среднее давление при этом падает, а максимальное падение происходит в центре трубы. Расход также испытывает осцилляции, но на гораздо меньших частотах и с незначительными амплитудами.

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

Настоящая работа развивает идеи, лежащие в основе [9]. Главным отличием является рассмотрение полной постановки задачи с учетом уравнения энергии и термодинамического уравнения состояния для совершенного газа, тогда как в [9] жидкость считалась слабо сжимаемой. Мотивация заключается в том, чтобы проследить изменения температуры в резонансных режимах. При численном моделировании также учитываются изменения вязкости и теплопроводности вместе с изменениями температуры. Однако эти изменения оказываются несущественными.

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

Рассмотрим течение вязкой сжимаемой теплопроводной жидкости (газа) в длинном плоском канале с жесткими стенками. Некоторые части стенок образуют вибрирующие секции, способные осциллировать с постоянной частотой $\omega $ в направлении, перпендикулярном основному потоку, рис. 1а. Амплитуда осцилляций $\varepsilon $ мала по сравнению с шириной канала $l$. Отдельная секция с соответствующими граничными условиями показана на рис. 1б.

Рис. 1.

(а) Геометрия плоского канала с пятью осциллирующими секциями стенок; (б) Отдельная осциллирующая секция.

Течение описывается двумерными уравнениями Навье–Стокса в декартовых координатах $\left( {x,y} \right)$, уравнением неразрывности и уравнением баланса энергии с коэффициентами $\mu $ динамической вязкости и ${{\theta }}$ теплопроводности, зависящими от температуры $T$. Система замыкается уравнениями состояния (термодинамическим и энергетическим) совершенного газа с постоянными теплоемкостями ${{c}_{p}}$ и ${{c}_{{\text{v}}}}$

(1.1)
$\frac{{\partial \rho }}{{\partial t}} + u\frac{{\partial \rho }}{{\partial x}} + \rho \frac{{\partial u}}{{\partial x}} + {v}\frac{{\partial \rho }}{{\partial y}} + \rho \frac{{\partial {v}}}{{\partial y}} = 0$
(1.2)
$\rho \frac{{Du}}{{Dt}} + \frac{{\partial p}}{{\partial x}} = \frac{\partial }{{\partial x}}\left[ {\mu \left( {2\frac{{\partial u}}{{\partial x}} - \frac{2}{3}\left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial {v}}}{{\partial y}}} \right)} \right)} \right] + \frac{\partial }{{\partial y}}\left[ {\mu \left( {\frac{{\partial u}}{{\partial y}} + \frac{{\partial {v}}}{{\partial x}}} \right)} \right]$
(1.3)
$\rho \frac{{D{v}}}{{Dt}} + \frac{{\partial p}}{{\partial y}} = \frac{\partial }{{\partial x}}\left[ {{{\mu }}\left( {\frac{{\partial u}}{{\partial y}} + \frac{{\partial {v}}}{{\partial x}}} \right)} \right] + \frac{\partial }{{\partial y}}\left[ {{{\mu }}\left( {2\frac{{\partial {v}}}{{\partial y}} - \frac{2}{3}\left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial {v}}}{{\partial y}}} \right)} \right)} \right]$
(1.4)
${{\rho }}\frac{{De}}{{Dt}} + p\,\left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial {v}}}{{\partial y}}} \right) = \frac{\partial }{{\partial x}}\left( {{{\theta }}\frac{{\partial T}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {{{\theta }}\frac{{\partial T}}{{\partial y}}} \right) + {{\mu }}\,\Theta $
(1.5)
$p = \rho RT,\quad e = {{c}_{{v}}}T$
(1.6)
$\mu = {{\mu }_{{{\text{contr}}}}}\frac{{{{T}_{{_{{{\text{contr}}}}}}} + {{T}_{\mu }}}}{{T + {{T}_{\mu }}}}{{\left( {\frac{T}{{{{T}_{{_{{{\text{contr}}}}}}}}}} \right)}^{{3{\text{/}}2}}},\quad \theta = {{\theta }_{{_{{{\text{contr}}}}}}}\frac{{{{T}_{{_{{{\text{contr}}}}}}} + {{T}_{\theta }}}}{{T + {{T}_{\theta }}}}{{\left( {\frac{T}{{{{T}_{{_{{{\text{contr}}}}}}}}}} \right)}^{{3{\text{/}}2}}}$

В системе (1.1)–(1.6) приняты следующие обозначения: $u,\;{v}$ – компоненты вектора скорости; $p$ – давление, ${{\rho }}$ – плотность, $T$ – абсолютная температура; $e$ – удельная внутренняя энергия, $R = {{c}_{p}} - {{c}_{{v}}}$ – газовая постоянная, а $\Theta $ – диссипативная функция

$\Theta = 2{{\left( {\frac{{\partial u}}{{\partial x}}} \right)}^{2}} + 2{{\left( {\frac{{\partial {v}}}{{\partial y}}} \right)}^{2}} + {{\left( {\frac{{\partial u}}{{\partial y}} + \frac{{\partial {v}}}{{\partial x}}} \right)}^{2}} - \frac{2}{3}{{\left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial {v}}}{{\partial y}}} \right)}^{2}}$

Уравнения (1.2) и (1.3) выписаны в рамках гипотезы Стокса – объемная вязкость считается равной нулю: ${{{{\zeta }}}^{{vol}}} = {{\lambda }} + 2{{\mu /}}3 \equiv 0$, где ${{\lambda }}$ – второй коэффициент вязкости.

Формулы (1.6) – зависимости Сазерленда коэффициентов вязкости и теплопроводности от температуры для совершенного газа, в которых ${{{{\mu }}}_{{_{{contr}}}}}$ и ${{{{\theta }}}_{{_{{contr}}}}}$ – контрольные вязкость и теплопроводность при контрольной температуре ${{T}_{{_{{contr}}}}}$, а ${{T}_{{{\mu }}}}$ и ${{T}_{\theta }}$ – постоянные Сазерленда, характеризующие конкретный газ.

Систему (1.1)–(1.6) необходимо дополнить начальными и граничными условиями. На входе в канал задаются постоянные входное давление и температура: $p = {{P}_{{in}}},$ $T = {{T}_{{in}}}$. На выходе из канала – постоянное выходное давление и “мягкое” условие по температуре: $p = {{P}_{{out}}},$ $\partial T{\text{/}}\partial x = 0$. На вибрирующих секциях стенок канала задаются гармонические осцилляции вертикальной скорости и условие теплоизоляции

(1.7а)
$\left\{ {x \in {{{\mathbf{X}}}_{{{\mathbf{vibe}}}}},y = 0 + \varepsilon \sin \left( {\omega t} \right)} \right\}:\quad u = 0,\quad {v} = {{V}_{d}}\left( t \right) = \varepsilon \omega \cos \left( {\omega t} \right),\quad \partial T{\text{/}}\partial y = 0$
(1.7б)
$\left\{ {x \in {{{\mathbf{X}}}_{{{\mathbf{vibe}}}}},y = l + \varepsilon \sin \left( {\omega t} \right)} \right\}:\quad u = 0,\quad {v} = {{V}_{u}}\left( t \right) = \varepsilon \omega \cos \left( {\omega t} \right),\quad \partial T{\text{/}}\partial y = 0$

На покоящихся частях стенок ставятся обычные условия прилипания и теплоизоляции. В условиях (1.7) ${{{\mathbf{X}}}_{{{\mathbf{vibe}}}}}$ означает множество тех $x \in \left[ {0;\Lambda } \right]$, которые соответствуют вибрирующим секциям. Здесь и далее, $\Lambda $ – длина всего канала, $L$ – длина одной секции.

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

(1.8)
$\left\{ {t = 0} \right\}:\quad u = 0,\quad {v} = 0,\quad p = {{p}_{0}},\quad T = {{T}_{0}}$

Замкнутая математическая модель (1.1)–(1.8) составляет основу настоящей статьи, состоящей из двух частей. В первой части аналитически рассматриваются линеаризованные решения для одной осциллирующей секции (рис. 1б). Во второй нелинейная система (1.1)–(1.8) решается конечно-разностным методом в полной геометрии канала с пятью осциллирующими секциями (рис. 1а).

2. ЛИНЕАРИЗОВАННЫЕ РЕШЕНИЯ

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

Из уравнений (1.1)–(1.3) можно исключить компоненты скоростей

(2.1)
$\frac{{{{\partial }^{2}}{{\rho }}{\kern 1pt} '}}{{\partial {{t}^{2}}}} = \frac{{a_{0}^{2}}}{{{\gamma }}}\Delta {{\rho }}{\kern 1pt} '\; + {{{{\rho }}}_{0}}R\Delta T{\kern 1pt} '\; + \frac{4}{3}{{\nu }_{0}}\frac{\partial }{{\partial t}}\Delta {{\rho }}{\kern 1pt} '$
где $a_{0}^{2} = {{\gamma }}R{{T}_{0}}$ – невозмущенная скорость звука, ${{\gamma }} = {{c}_{p}}{\text{/}}{{c}_{{v}}}$ – показатель адиабаты, ${{\nu }_{0}}$ – невозмущенная кинематическая вязкость.

Тогда уравнение (2.1) и линеаризованное уравнение (1.4) образуют замкнутую систему для определения отклонений плотности и температуры. Отклонение давления определяется первым (линеаризованным) уравнением (1.5) в явном виде.

Исключая из уравнений (2.1) и (1.4) температуру, можно прийти к следующему уравнению для плотности (штрихи в дальнейшем изложении опустим)

(2.2)
$\frac{\partial }{{\partial t}}\left[ {\frac{{{{\partial }^{2}}}}{{\partial {{t}^{2}}}} - a_{0}^{2}\Delta } \right]{{\rho }}{\kern 1pt} '\; - \Delta \left[ {\left( {b_{0}^{2} + {{\gamma }}c_{0}^{2}} \right)\frac{{{{\partial }^{2}}}}{{\partial {{t}^{2}}}} - a_{0}^{2}c_{0}^{2}\Delta - {{\gamma }}b_{0}^{2}c_{0}^{2}\frac{\partial }{{\partial t}}\Delta } \right]{{\rho }}{\kern 1pt} ' = 0$
где $b_{0}^{2} = 4{{\nu }_{0}}{\text{/}}3$, а $c_{0}^{2} = {{{{\theta }}}_{0}}{\text{/}}{{{{\rho }}}_{0}}{{c}_{{\text{p}}}}$ – температуропроводность жидкости.

Отметим, что при отсутствии диссипативных эффектов (вязкости и теплопроводности) уравнение (2.2) сводится к обычному волновому уравнению. Поскольку общее решение уравнения (2.2) достаточно громоздко, более наглядным представляется разделить исследование влияния диссипативных эффектов: сначала рассмотреть вязкую, но нетеплопроводную жидкость ($c_{0}^{2} = 0$), а затем невязкую, но теплопроводную. Оба результата сравнить затем со случаем идеальной (невязкой и нетеплопроводной) жидкости. Можно показать, что температура в обоих случаях удовлетворяет уравнению, полностью аналогичному уравнению (2.2) для плотности.

Граничные условия для уравнения (2.2) примем в следующем виде

(2.3а)
$x = 0:\quad {{\rho }} = {{{{\rho }}}_{{in}}},\quad x = L:\quad \rho = {{\rho }_{{\text{0}}}}$
(2.3б)
$y = 0:\quad {{\rho }} = {{{{\rho }}}_{0}}\left[ {1 + \left( {S{\text{/}}{{{{\rho }}}_{0}}} \right){\text{sin}}\left( {\omega t} \right)} \right],\quad y = l:\quad {{\rho }} = {{{{\rho }}}_{0}}\left[ {1 + \left( {S{\text{/}}{{{{\rho }}}_{0}}} \right){\text{sin}}\left( {\omega t} \right)} \right]$
где $S$ – некоторая константа.

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

Начальные условия имеют вид

(2.4)
$t = 0:\quad {{\rho }} = {{{{\rho }}}_{0}},\quad \partial {{\rho /}}\partial t = 0$

Для решения уравнения (2.2) применим методом Фурье разделения переменных. Изложим только общую схему без вычислительных подробностей.

1. Введением переменной ${{\rho }} = {{{{\rho }}}_{0}} + S\sin \left( {\omega \,t} \right) + {{\tilde {\rho }}}$ обнуляются граничные условия (2.3б) в направлении y. Переменные разделяются. Решение раскладывается по собственным функциям $\sin \left( {\pi ny{\text{/}}l} \right)$: ${{\tilde {\rho }}} = \sum\limits_{n = 0}^\infty {{{F}_{n}}\left( {x,t} \right) \cdot \sin \frac{{\pi ny}}{l}} $

2. Введением новой переменной ${{F}_{n}} = {{\tilde {F}}_{n}} - \left( {{{{{\rho }}}_{{in}}} - {{{{\rho }}}_{0}}} \right){{c}_{n}}\left( {x - L} \right){\text{/}}L - S{{c}_{n}}\sin \left( {\omega t} \right)$ обнуляются граничные условия (2.3а) в направлении $x$. Переменные разделяются. Решение раскладывается по собственным функциям $\sin \left( {\pi kx{\text{/}}L} \right)$: ${{\tilde {F}}_{n}} = \sum\limits_{k = 0}^\infty {{{T}_{{nk}}}\left( t \right) \cdot \sin \frac{{\pi kx}}{L}} $.

3. Задача сводится к обыкновенному дифференциальному уравнению для неизвестных функций времени ${{T}_{{nk}}}\left( t \right)$ с заданными начальными условиями.

Приведем только окончательные уравнения. При отсутствии теплопроводности функции ${{T}_{{nk}}}\left( t \right)$ должны удовлетворять следующему обыкновенному дифференциальному уравнению

(2.5а)
$\frac{{{{d}^{2}}{{T}_{{nk}}}}}{{d{{t}^{2}}}} + 2{{\alpha }_{{nk}}}\frac{{d{{T}_{{nk}}}}}{{dt}} + \omega _{{0nk}}^{2}{{T}_{{nk}}} = A_{{nk}}^{'}\sin \left( {\omega \,t} \right) + B_{{nk}}^{'}\cos \left( {\omega t} \right) + C_{{nk}}^{'}$
$2{{\alpha }_{{nk}}} = b_{0}^{2}{{\Omega }_{{nk}}},\quad \omega _{{0nk}}^{2} = a_{0}^{2}{{\Omega }_{{nk}}},\quad {{\Omega }_{{nk}}} = {{\left( {\frac{{\pi n}}{l}} \right)}^{2}} + {{\left( {\frac{{\pi k}}{L}} \right)}^{2}}$
$A_{{nk}}^{'} = S{{\left( {\frac{{\pi n}}{l}} \right)}^{2}}{{c}_{n}}{{d}_{k}}\frac{{a_{0}^{2}}}{{{\gamma }}},\quad B_{{nk}}^{'} = S{{\left( {\frac{{\pi n}}{l}} \right)}^{2}}{{c}_{n}}{{d}_{k}}b_{0}^{2}\omega ,\quad C_{{nk}}^{'} = \frac{{a_{0}^{2}}}{{{\gamma }}}{{\left( {\frac{{\pi n}}{l}} \right)}^{2}}{{c}_{n}}{{e}_{k}}\left( {{{{{\rho }}}_{{in}}} - {{{{\rho }}}_{0}}} \right)$

При отсутствии вязкости уравнение для функций ${{T}_{{nk}}}\left( t \right)$ принимает вид

(2.5б)
$\frac{{{{d}^{3}}{{T}_{{nk}}}}}{{d{{t}^{3}}}} + {{\gamma }}{{\beta }_{{nk}}}\frac{{{{d}^{2}}{{T}_{{nk}}}}}{{d{{t}^{2}}}} + \omega _{{0nk}}^{2}\frac{{d{{T}_{{nk}}}}}{{dt}} + {{{{\beta }}}_{{nk}}}\omega _{{0nk}}^{2}{{T}_{{nk}}} = A_{{nk}}^{{''}}\sin \left( {\omega t} \right) + B_{{nk}}^{{''}}\cos \left( {\omega \,t} \right) + C_{{nk}}^{{''}}$
${{{{\beta }}}_{{nk}}} = c_{0}^{2}{{\Omega }_{{nk}}},\quad \omega _{{0nk}}^{2} = a_{0}^{2}{{\Omega }_{{nk}}},\quad A_{{nk}}^{{''}} = S{{\left( {\frac{{\pi n}}{l}} \right)}^{2}}{{c}_{n}}{{d}_{k}}c_{0}^{2}\left( {a_{0}^{2}{{{\left( {\frac{{\pi n}}{l}} \right)}}^{2}} - {{\gamma }}{{\omega }^{2}}} \right)$
$B_{{nk}}^{{''}} = S{{\left( {\frac{{\pi n}}{l}} \right)}^{2}}{{c}_{n}}{{d}_{k}}a_{0}^{2}\omega ,\quad C_{{nk}}^{{''}} = a_{0}^{2}c_{0}^{2}{{\left( {\frac{{\pi n}}{l}} \right)}^{4}}{{c}_{n}}{{e}_{k}}\left( {{{{{\rho }}}_{{in}}} - {{{{\rho }}}_{0}}} \right)$

Коэффициенты Фурье ${{c}_{n}}$, ${{d}_{k}}$ и ek в уравнениях (2.5) определяются соотношениями

(2.6)
${{c}_{n}} = \frac{2}{{\pi n}}\left( {1 - {{{\left( { - 1} \right)}}^{n}}} \right),\quad {{d}_{k}} = \frac{2}{{\pi k}}\left( {1 - {{{\left( { - 1} \right)}}^{k}}} \right),\quad {{e}_{k}} = \frac{2}{{\pi k}}$

Начальные условия для уравнений (2.5) имеют вид:

(2.7)
$t = 0:\quad {{T}_{{nk}}} = \left( {{{{{\rho }}}_{{in}}} - {{{{\rho }}}_{0}}} \right){{c}_{n}}{{e}_{k}},\quad \frac{{d{{T}_{{nk}}}}}{{dt}} = 0$

Решение обыкновенного дифференциального уравнения (2.5а), отвечающего случаю отсутствия теплопроводности, является суммой общего (2.8а) и частного (2.9а) решений

(2.8а)
$T_{{nk,vis}}^{{gen}}\left( t \right) = {{e}^{{ - {{\alpha }_{{nk}}}t}}}\left[ {D_{{1nk}}^{'}\sin (\sqrt {\omega _{{0nk}}^{2} - {{\alpha }}_{{nk}}^{2}} \cdot t) + D_{{2nk}}^{'}\cos (\sqrt {\omega _{{0nk}}^{2} - {{\alpha }}_{{nk}}^{2}} \cdot t)} \right]$
(2.9а)
$T_{{nk,vis}}^{{part}}\left( t \right) = \Psi _{{nk}}^{'}\sin (\omega t + {{\theta }}_{{nk}}^{'}) + \frac{{C_{{nk}}^{'}}}{{\omega _{{0nk}}^{2}}}$

Аналогично, решение обыкновенного дифференциального уравнения (2.5б), отвечающего случаю отсутствия вязкости, является суммой общего (2.8б) и частного (2.9б) решений, однако найти собственные значения аналитически в этом случае не получается, так как они удовлетворяют алгебраическому уравнению третьей степени. При этом легко можно показать, что действительные части всех трех собственных чисел отрицательны. Таким образом, неограниченного растущих во времени решений у уравнения (2.5б) нет

(2.8б)
$T_{{{\text{nk,cond}}}}^{{{\text{gen}}}}\left( t \right) = D_{{1nk}}^{{''}}{{e}^{{{{\lambda }_{1}}t}}} + D_{{2nk}}^{{''}}{{e}^{{{{\lambda }_{2}}t}}} + D_{{3nk}}^{{''}}{{e}^{{{{\lambda }_{3}}t}}},\quad {{\lambda }_{i}} = {{\left( {{{\lambda }_{x}}} \right)}_{i}} + i{{\left( {{{\lambda }_{y}}} \right)}_{i}},\quad {{\left( {{{\lambda }_{x}}} \right)}_{i}} < 0,\quad i = 1,2,3$
(2.9б)
$T_{{nk,cond}}^{{part}}\left( t \right) = \Psi _{{nk}}^{{''}}\sin (\omega t + {{\theta }}_{{nk}}^{{''}}) + \frac{{C_{{nk}}^{{''}}}}{{{{{{\beta }}}_{{nk}}}\omega _{{0nk}}^{2}}}$
$\Psi _{{nk}}^{{''}} = \frac{{\sqrt {A_{{nk}}^{{''2}} + B_{{nk}}^{{''2}}} }}{{\sqrt {\beta _{{nk}}^{2}{{{(\omega _{{0nk}}^{2} - \gamma {{\omega }^{2}})}}^{2}} + {{\omega }^{2}}{{{(\omega _{{0nk}}^{2} - {{\omega }^{2}})}}^{2}}} }},\quad {\text{tg}}{{\theta }_{{nk}}} = \frac{{{{B}_{{nk}}}{{\beta }_{{nk}}}(\omega _{{0nk}}^{2} - \gamma {{\omega }^{2}}) - {{A}_{{nk}}}\omega (\omega _{{0nk}}^{2} - {{\omega }^{2}})}}{{{{A}_{{nk}}}{{\beta }_{{nk}}}(\omega _{{0nk}}^{2} - \gamma {{\omega }^{2}}) + {{B}_{{nk}}}\omega (\omega _{{0nk}}^{2} - {{\omega }^{2}})}}$

Неизвестные константы $D_{{ink}}^{'},\;i = 1,2$ в общем решении (2.8а) и $D_{{ink}}^{{''}},\;i = 1,2,3$ в решении (2.8б) в принципе могут быть найдены из начальных условий (2.7). После этого выражения (2.8) и (2.9) дадут окончательное решение исходной начально-краевой задачи (2.2)–(2.4). Интереснее другое, определение условий резонанса – условий, при которых амплитуды частных $\Psi _{{nk}}^{'} = \Psi _{{nk}}^{{vis}}$ и $\Psi _{{nk}}^{{''}} = \Psi _{{nk}}^{{cond}}$ решений (2.9) достигают максимума.

В случае вязкой, но нетеплопроводной жидкости частоту резонанса удается получить аналитически. Дифференцированием функции $\Psi _{{nk}}^{{vis}}({{\omega }^{2}})$ можно показать, что наибольшая амплитуда частного решения (2.9а) достигается при

(2.10)
$\omega _{{res}}^{{vis}} = \frac{{3a_{0}^{2}}}{{4{{\nu }_{0}}}}\sqrt {\sqrt {1 + \frac{{32\nu _{0}^{2}}}{{9a_{0}^{2}}}\left( {{{{\left( {\frac{{\pi n}}{l}} \right)}}^{2}} + {{{\left( {\frac{{\pi k}}{L}} \right)}}^{2}}} \right)} - 1} $

В случае невязкой, но теплопроводной жидкости дифференцированием функции $\Psi _{{nk}}^{{cond}}({{{{\omega }}}^{2}})$ можно показать, что квадрат частоты резонанса ${{({{\omega }}_{{res}}^{{cond}})}^{2}}$ удовлетворяет алгебраическому уравнению четвертой степени (вследствие громоздкости, уравнение не приводится). Поиск корней этого уравнения осуществлялся численно.

В случае отсутствия вязкости и теплопроводности (идеальная жидкость) выражение для частоты резонанса заметно упрощается

(2.11)
${{\omega }}_{{res}}^{{ideal}} = {{\omega }}_{{nk}}^{0} = {{a}_{0}}\sqrt {{{{\left( {\frac{{\pi n}}{l}} \right)}}^{2}} + {{{\left( {\frac{{\pi k}}{L}} \right)}}^{2}}} $

Однако поведение частного решения при этом существенно меняется – оно начинает неограниченно расти со временем:

$T_{{nk,ideal}}^{{part}}\left( t \right) = - \frac{{A_{{nk}}^{'}}}{{2{{{{\omega }}}_{{0nk}}}}}t \cdot \cos \left( {{{{{\omega }}}_{{0nk}}}t} \right) + \frac{{C_{{nk}}^{'}}}{{{{\omega }}_{{0nk}}^{2}}}$

Анализ полученных решений проведем для воздуха. Длину секции возьмем равной длине канала. Наиболее важным представляется вывод, что величины “вязкой” (2.10) и “теплопроводной” частот резонанса отличаются от соответствующей частоты (2.11) для идеальной жидкости менее, чем на сотою долю процента. Однако, несмотря на практическое совпадение всех трех частот, в поведении соответствующих решений есть существенное отличие: наибольшая амплитуда в случае идеальной жидкости стремится к бесконечности. Наличие же диссипативных эффектов (вязкости или теплопроводности) приводит к тому, что наибольшая амплитуда сразу становится конечной. Анализ показывает, что наиболее опасной (ведущей) частотой является частота ${{({{\omega }}_{{res}}^{{lead}})}_{{11}}}$, соответствующая первой гармонике $n = k = 1$. Для более высоких гармоник усиление амплитуды падает экспоненциально. Таким образом, резонанс возможен только на нескольких первых частотах.

Зависимость общего решения (плотности) от продольной координаты $x$ при $y = 0.05\,{\text{м}}$ в момент $t = 0.2\,{\text{c}}$ в вязком нетеплопроводном случае при различных частотах из резонансного спектра показано на рис. 2. Видно, что наибольшей амплитуды плотность достигает на ведущей частоте ${{\omega }}_{{11}}^{{res}}$. С другой стороны, решение для четной частоты ${{\omega }}_{{22}}^{{res}}$ испытывает минимальные осцилляции, незаметные в выбранном масштабе.

Рис. 2.

Зависимость плотности от координаты x в вязком нетеплопроводном случае при $y = 0.05$ м в момент $t = 0.2$ c для различных частот из резонансного спектра.

3. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ

Перейдем к численному решению исходной нелинейной системы уравнений (1.1)–(1.8) движения вязкого сжимаемого теплопроводного газа в полной геометрии канала с пятью осциллирующими секциями (рис. 1а). Основным методом интегрирования будет метод разделения по физическим процессам [10]. Интегрирование по времени разбивается на два этапа: на первом (невязком) этапе уравнения (1.1)(1.6) решаются при отсутствии вязкости и теплопроводности, на втором – учитывается действие этих диссипативных слагаемых на параметры, рассчитанные на первом этапе.

При отсутствии вязкости и теплопроводности система (1.1)–(1.6) сводится к уравнениям газовой динамики. Для их решения применялся другой вид метода расщепления – разделение по направлениям – первый этап интегрирования также разбивался на два промежуточных подшага. На первом подшаге в уравнениях газовой динамики остаются слагаемые, содержащие производные только по координате x, а на втором – только по координате y. Получившиеся уравнения решались по явным схемам первого порядка точности, основанным на представлении уравнений в характеристической форме и исключающим возможные численные осцилляции.

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

Расчетная область (прямоугольный канал) $\Omega = \left[ {0;\Lambda } \right] \times \left[ {0;l} \right]$ разбивается на прямоугольные ячейки: четыреста ячеек в направлении x (${{N}_{x}} = 400$) и двести в направлении y (${{N}_{y}} = 200$). Изначально канал заполнен покоящейся жидкостью, давление и температура постоянны. В начальный момент t = 0 давление на входе поднимается, одновременно начинаются осцилляции стенок.

Важно отметить, что стенки во всех расчетах покоятся, а их влияние на течение идет через граничные условия (1.7а) и (1.7б). В качестве обоснования подобного подхода можно сказать, что реальное движение стенок (осцилляции) покрывает меньшее расстояние, чем ширина одной ячейки.

В качестве тестовых задач рассматривались стационарное течение Куэтта вязкой сжимаемой теплопроводной жидкости и нестационарное течения Пуазейля несжимаемой и нетеплопроводной жидкости. В обычных предположениях обе задачи допускают аналитические решения. Для всех задач получено хорошее соответствие численных и аналитических результатов – ошибка составляла менее одного процента. Также проводилось исследование на сходимость по сетке. Обе тестовые задачи рассчитывались на грубой (${{N}_{x}} = 300,{{N}_{y}} = 175$) и сгущенной сетках (${{N}_{x}} = 500$, ${{N}_{y}} = 225$). Распределения основных параметров (скорости, температуры) отличались менее чем на один процент. Условие резонанса (2.11) оказывалось нечувствительно к конкретному виду сетки.

Перейдем к результатам расчетов течения в канале с осциллирующими секциями. Все пять секций имеют одинаковую длину $L = 0.1$ м, длина канала $\Lambda = 1$ м, ширина $l = 0.1$ м. В качестве рабочей жидкости (газа) используем газ в пять раз более вязкий, чем воздух. Установление течения в рассматриваемых условиях происходит примерно за 60 с. Важность выхода на стационарный режим заключается в следующем: мы должны четко показать, что кумулятивный резонансный эффект (рост) расхода в установившемся режиме при наличии осциллирующих секций существенно превосходит установившийся расход в канале с покоящимися стенками.

Ключевыми безразмерными параметрами расчетов будут вибрационное число Рейнольдса и число Маха. Числа Маха будут крайне малы M ~ 0.1, однако вибрационные числа Рейнольдса, соответствующие резонансным режимам, оказываются достаточно большими $\alpha \sim 1.5 \times {{10}^{5}}$. Основным размерным параметром, описывающим течение, будет массовый расход рабочего газа Q. В качестве контрольного сечения выбиралось выходное сечение канала. Давление на выходе канала всегда остается атмосферным.

Отметим, что упрощенная аналитическая постановка разд. 2, в которой граничные условия (2.3) первого рода (типа Дирихле) только для плотности однородны, приводит к тому, что наступление резонанса возможно на ведущей частоте $n = 1,$ $k = 1$. Однако полная нелинейная постановка включает неоднородные граничные условия первого и второго рода (типа Неймана) по разным параметрам (давлению, продольной скорости и температуре). В результате наступление резонанса в расчетах становится возможным и на более низких частотах: $n = 1,$ $k = 0$ и $n = 0,$ $k = 1$.

Во всех проведенных расчетах резонанс наступал при частоте, соответствующей $n = 0,$ $k = 1$, потому что она оказывалось самой низкой – ${{\omega }}_{{01}}^{{res}} = 1089$ Гц. Заметим, что это значение представляется достаточно реалистичным для технических приложений. Усиление амплитуды на этой частоте (назовем ее ведущей) оказывалось наиболее интенсивным. Кроме того, кумулятивный резонанс (по истечении некоторого времени) наблюдался только на ведущей частоте ${{\omega }}_{{01}}^{{res}}$.

При наличии осциллирующих секций все параметры течения начинают периодически меняться вместе с осцилляциями. Однако, чем больше разница между частотой осцилляций стенок и частотой резонанса, тем меньше амплитуда этих колебаний параметров. На рис. 3 и рис. 4 показано поведение давления и температуры в обычном и резонансном режимах в момент $t = 40$ с. Верхние и нижние секции осциллируют на одной частоте в синфазном режиме. Максимальная амплитуда осцилляций поперечной скорости составила ${{V}_{0}} = 0.05$ м/с или примерно ${{10}^{{ - 2}}} \cdot {{U}_{0}}$.

Рис. 3.

Распределение давления в канале с пятью осциллирующими секциями при $t = 40\,{\text{c}}$: (а) осцилляции стенок не приводят к кумулятивному резонансу, частота ${{\omega }}_{{11}}^{{res}}$; (б) осцилляции стенок дают кумулятивный резонанс на ведущей частоте ${{\omega }}_{{01}}^{{res}}$.

Рис. 4.

Распределение температуры в канале с пятью осциллирующими секциями при $t = 40\,{\text{c}}$: (а) осцилляции стенок не приводят к кумулятивному резонансу, частота ${{\omega }}_{{11}}^{{res}}$; (б) осцилляции стенок дают кумулятивный резонанс на ведущей частоте ${{\omega }}_{{01}}^{{res}}$.

На рис. 3а и рис. 4а осцилляции стенок происходят на частоте из резонансного спектра $\omega _{{{\text{11}}}}^{{{\text{res}}}}$, не приводящей к кумулятивному эффекту. Хорошо заметен периодический характер решения. На рис. 3б и рис. 4б стенки осциллируют на ведущей частоте $\omega _{{{\text{01}}}}^{{{\text{res}}}}$. Виден резкий рост давления по центру канала, соответствующий результатам линейной теории (рис. 2). Жидкость разогревается к середине канала, однако изменения температуры составляют не более градуса.

Главный результат, демонстрирующий многократное возрастание массового расхода, показан на рис. 5. Только на ведущей резонансной частоте $\omega _{{{\text{01}}}}^{{{\text{res}}}}$ после некоторых “установившихся” периодических изменений давления, плотности и температуры возникает кумулятивный резонанс – резкое увеличение суммарного массового расхода (черная кривая). Видно, что кумулятивный резонанс происходит только на ведущей частоте примерно через минуту после начала осцилляций. Расход для канала с покоящимися стенками к этому моменту уже вышел на установившийся режим. В околорезонансном случае ${{\omega }} = 1.1{{\omega }}_{{01}}^{{res}}$ эффект уже значительно меньше или вовсе отсутствует ${{\omega }} = 0.9{{\omega }}_{{01}}^{{res}}$. Для более высоких частот резонансного спектра ${{\omega }}_{{10}}^{{res}}$ и ${{\omega }}_{{11}}^{{res}}$ расчеты показывают, что разница по сравнению с отсутствием осцилляций практически не наблюдается.

Рис. 5.

Зависимость среднего массового расхода на выходе канала от времени: 1 – покоящиеся стенки; 2 – осцилляции секций на ведущей резонансной частоте; 3, 4 – осцилляции на “околорезонансных” частотах.

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

Отметим, что изменения температуры даже в резонансном режиме являются ничтожными (полградуса), что не приводит к существенному изменению вязкости и теплопроводности, которые остаются практически постоянными. Таким образом, кумулятивный эффект нарастания расхода является следствием именно резонансного эффекта.

ЗАКЛЮЧЕНИЕ

Основным результатом работы является обнаруженный на основе линеаризованных решений и численного моделирования кумулятивный резонанс – значительное увеличение массового расхода жидкости при заданном постоянном перепаде давления. Получено точное аналитическое выражение для частоты резонанса в канале с осциллирующими стенками в случае идеальной (невязкой и нетеплопроводной) жидкости. Учет вязкости и теплопроводности приводит к ее незначительному (менее чем на сотые доли процента) изменению. Однако амплитуда решения сразу становится конечной.

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

Основное, но далеко не единственное приложение полученных результатов – течение газа и нефти в трубопроводах или воды в помповых насосах.

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

  1. Ганиев О.Р., Ганиев Р.Ф., Звягин А.В., Украинский Л.Е., Устенко И.Г. Повышение нефтеотдачи пластов на основе волноводных эффектов. Немонотонность затухания двумерных волн в волноводе // Спр. инж. журн. 2016. 228. № 3. С. 42–48.

  2. Ганиев О.Р., Ганиев Р.Ф., Звягин А.В., Украинский Л.Е., Устенко И.Г. Повышение нефтеотдачи пластов на основе волноводных эффектов. Распределение в волноводе силы, действующей на частицы в гармоническом волновом поле // Спр. инж. журн. 2016. 229. № 4. С. 49–56.

  3. Secomb T.W. Flow in a channel with pulsating walls // J. Fluid. Mech. 1978. V. 88. P. 273–288.

  4. Hall P., Papageorgiou D.T. The onset of chaos in a class of Navier – Stokes solutions // J. Fluid. Mech. 1999. V. 393. P. 59–87.

  5. Mingalev S.V., Filippov L.O., Lyubimova T.P. Flow rate in a channel with small-amplitude pulsating walls // Euro. J. Mech. B/Fl. 2015. V. 51. P. 1–7.

  6. Тукмаков А.Л. Течение газа в плоском канале с вибрирующими стенками // Инж.-физ. журн. 2002. Т. 75. № 6. С. 109–115.

  7. Smirmov N.N., Shugan I.V., Legros J.C. Streaming flows in a channel with elastic walls // Phys. Fluids. 2002. V. 14. № 10. P. 1–10.

  8. Bagchi S., Gupta K., Kushari A., Iyengar N.G.R. Experimental study of pressure fluctuations and flow perturbations in air flow through vibrating pipes // J. Sound. Vibr. 2009. V. 328. P. 441–455.

  9. Logvinov O.A., Malashin A.A. Resonance phenomena in a channel with oscillating walls // Euro. J. Mech. B/Fl. 2020. 83. 205–211.

  10. Ковеня В.М., Яненко Н.Н. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, 1981. 304 с.

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