CC BY
43
Вычислительные технологии
Том 13, № 5, 2008
Численное моделирование волн Кельвина—Гельмгольца в слабо неравновесном молекулярном газе*
Ю. Н. Григорьев Институт вычислительных технологий СО РАН, Новосибирск, Россия
email: grigor@ict.nsc.ru
И. В. Ершов
Новосибирская государственная академия водного транспорта, Россия
email: i_ershov@ngs.ru
К. И. Зырянов
Новосибирский государственный архитектурно-строительный университет,
Россия
An influence of the bulk viscosity on suppression of the Kelvin—Helmholtz instability in time-dependent shear layers was investigated in the frameworks of Navier—Stocks equations of compressible gas.
Введение
В работах [1-3] показано, что присутствие термического возбуждения заметно усиливает подавление вихревых возмущений в молекулярном газе. Рассматривалась эволюция вихревого возмущения в виде вихря Рэнкина, наложенного на сдвиговое течение Куэтта. Несмотря на простоту, такая модель при соответствующем выборе параметров вполне удовлетворительно воспроизводит процесс диссипации крупных вихревых структур, возникающих в нелинейной стадии ламинарно-турбулентного перехода и в развитых турбулентных течениях. Вместе с тем, как отмечалось уже в [1], модель не охватывает начальной стадии генерации структур малыми возмущениями, усиливающимися за счет кинетической энергии несущего (среднего) потока. Это не позволяет более полно оценить влияние термической неравновесности на весь процесс эволюции вихревых возмущений.
C этой целью в данной работе в рамках полных уравнений Навье—Стокса сжимаемого теплопроводного газа рассматривается нелинейное развитие неустойчивости Кельвина—Гельмгольца вплоть до формирования и последующей диссипации структуры cat's-eye в эволюционирующем во времени слое сдвига с точкой перегиба на профиле скорости. Термическая неравновесность внутренних степеней свободы молекул газа варьировалась изменением отношения коэффициентов сдвиговой и объемной вязкости в тензоре напряжений. В качестве начальных возмущений использовались наиболее
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-01-00359).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.
растущие моды, рассчитанные в линейном приближении уравнений невязкой газовой динамики.
1. Постановка задачи
1.1. Основные уравнения и начально-краевые условия
В координатной плоскости х,у рассматривается сдвиговое течение, в котором основной (несущий) поток равномерной плотности д0 и температуры Т0 направлен вдоль оси абсцисс х и имеет профиль скорости
ив(у) = и
где 80 — длина, на которой наклон профиля скорости максимальный.
В момент времени Ь = 0 на основной поток накладывается двумерное возмущение с длиной волны Л и волновым вектором к = (в, 0). Динамика возмущенного течения описывается системой полных уравнений Навье—Стокса сжимаемого вязкого теплопроводного газа. Характерная длина для обезразмеривания определялась максимальным наклоном профиля скорости
¿0 = ио
'¿и3
¿у у=о_
-1
В качестве других масштабирующих величин были выбраны скорость и0, плотность д0 и температура Т0 и образованные из них характерные время т0 = 50/и0 и давление Р0 = в0 и^О. В определенных на их основе безразмерных переменных система уравнений Навье—Стокса имеет вид
¿в
+ д
дих + диу
дх д2и
дУ
0,
¿их др 1 / д2их
в ¿Ь дх Бе \ дх2 ду2
¿иу др 1 / д2иу 1 д 2и
¿Ь ду Бе \ дх2
в-
+
+
в+ 7(7 - 1)М2р
ду2
(1 + За) д / дих ЗБе дх\ дх
+ (1 + За) ^/'дих + ду дх
диу
дУ
диу
7(7 - 1)М2
Бе
дих + диу
+
ду дх
7 М2р = дТ.
дх (1
д
ЗБе
7
БеРг 2
+
д
д2Т ^ +
дх2
ду
д2Т \ ду2 )
дих
(1)
+
диу
дх ду
д
1Ь дЬ + Чх дх + Чу ду'
Критерии, входящие в систему (1), определены следующим образом. Число Рей-нольдса есть Бе = и080д0число Маха несущего потока равно М = иО/\/7^Т0, Рг = ^ср/Л0 — число Прандтля. Здесь Я — газовая постоянная, 7 = ср/си — показатель адиабаты, ср и с — удельные теплоемкости при постоянных давлении и объеме
1
2
2
соответственно, Ао — коэффициент теплопроводности. Коэффициент а равен отношению объемной вязкости к сдвиговой а = ^и характеризует степень неравновесности внутренних степеней свободы молекул газа. Предполагалось, что диссипативные коэффициенты в системе (1) не зависят от температуры и постоянны.
Обезразмеренный несущий поток задается соотношениями
и(у)^апЬ у, т = д. = 1, р. = -М. (2)
Краевая задача первоначально ставится в бесконечной полосе, центр которой совпадает с началом координат: х £ [—х0,х0]; у £ (—го, го). Ширина полосы по координате х выбирается равной длине волны возмущения А = 2п/в. При этом х0 = п/в. Асимптотические условия при у —> переносятся на у = ±у0, где ордината у0 определяется из условия 11апЬ у0 — 1| < 10-9. В расчетах принято у0 = 10.
При г > 0 на границах х = ±х0, у £ [—у0, у0] ставились периодические условия
их (г, х0, у) = — х0, у), Пу (£, х0, у) = Пу (£, — х0, у), е(г, х0, у) = —х0, у), р(г, х0, у) = р(г, — х0, у).
Соответственно, при у = ±у0, х £ [—х0, х0] принимались условия невозмущенного потока:
их (г, х, у0) = —их (г, х, — у0) = 1, Пу (г, х, у0) = Пу (г, х, —у0) = 0,
^(г, х, у0) = ^(г, х, — у0) = 1, р(г, х, у0) = р(г, х, —у0) = 1/7М2.
Начальные условия для поля скорости и термодинамических величин, включающие возмущения, определялись в виде
их(0, х, у) = 1апЬу + иХ(0, х, у), иу(0, х, у) = и'у(0, х, у), (3)
£?(0, х, у) = 1 + ^'(0, х, у), Т(0, х, у) = 1 + Т'(0, х, у), р(0, х, у) = —+ р'(0, х, у).
7М2
Использованные в расчетах возмущения компонент вектора скорости п'х, п'у и термодинамических величин Т', р' рассчитывались предварительно на основе результатов работы [4].
1.2. Расчет начальных возмущений
В качестве накладываемых на основной поток начальных возмущений выбраны параметры собственных линейных невязких колебаний с наибольшими коэффициентами нарастания, исследованные в [4]. Они определяются как решения системы уравнений газовой динамики [5], линеаризованной относительно стационарного основного течения (2):
д«х тт дих ' ди. др'
+ и.+ и
дг дх у ду дх'
диУ + ^ диУ = — дР, (4)
дг 5 дх ду'
м2^ дР+и. дР 1+ди^+диу -0
дг дх / дх ду Р = (е' + Т ')/7 М2.
Таблица 1. Волновые числа и фазовые скорости наиболее растущих невязких возмущений
м в С. в с.
0.0 0.445 0.427 0.190
0.2 0.426 0.426 0.181
0.5 0.397 0.356 0.141
0.7 0.326 0.309 0.101
Искомые решения представляются в виде плоских волн, распространяющихся в направлении основного потока
«, иУ, д', ТV} = {«(у), у(у), д(у), Т(у),р}ехр [г^(ж - сЪ)], (5)
где и, V, д, Т, р — амплитудные функции возмущений, г — мнимая единица, с = сг + гс. — комплексная фазовая скорость волны.
Подстановка (5) в систему (4) приводит к спектральной задаче, которая в [4] при сг = 0 была преобразована к скалярному уравнению для р :
сРр 2 весЬ у ¿р 2
--в
¿у2 ^апЬ у — гс.) ¿у
1 — М2( 1апЬу — гс.)2 р = 0, (6)
р I = р I , = 0.
1 1у = — те 1 1у = + те
Здесь р = рг + гр., рг(у), р.(у) — действительная и мнимая части комплексной собственной функции давления, собственными значениями задачи (6) являются волновые числа в.
Исходя из (6) в работе [4] для дозвукового несущего потока при М < 1 были найдены волновые числа в, фазовые скорости с и коэффициенты усиления вс наиболее растущих невязких мод. Их значения для различных чисел Маха, использованные в наших расчетах, приведены в табл. 1.
Кроме спектральных параметров для постановки начальных условий (3) необходимы амплитудные функции возмущений (5) — собственные функции спектральной задачи, которые в данной работе рассчитывались путем численного интегрирования задачи (6) для значений параметров из табл. 1. Уравнение (6) разделялось на действительную и мнимую части, которые приводились к нормальной системе уравнений первого порядка. Данные Коши при у = 0 для нее выбирались аналогично работе [4]. Система интегрировалась на промежутке [—уо,уо] с использованием процедуры Рунге—Кутты четвертого порядка. Другие амплитудные функции определялись через найденные функции р с помощью следующих из (4), (5) соотношений
г (1апЬ у + %Сг \ { 2 . Д г (1апЬ у + гсА ¿р
и = — -о-о V весь у + гвр , V = — -к-- —,
в \tanh2 у + с2г)\ У ' в \tanh2 у + с2) ¿у'
в = М2р , Т = (7 — 1)М2р .
В результате вводимые в расчет начальные возмущения (3) определялись как вещественная часть рассчитанных решений (5):
я'(ж, у, Ъ) = И,е 4(у)ехр(¿вж) ехр (вс.Ъ),
где q' = (и^, и'у, д', Т', р'), я(у) = (и, V, в , Т, р). Следует заметить, что большим волновым числам (меньшим числам Маха) соответствуют большие амплитуды вводимых возмущений.
2. Численные расчеты эволюции возмущений
2.1. Разностная схема
Расчеты велись на регулярной сетке с одинаковым шагом к по обеим пространственным переменным. Система (1) аппроксимировалась весовой конечно-разностной схемой с расщеплением по физическим процессам и пространственным переменным [6], которая уже использовалась в предыдущих работах [1-3]. В операторной форме схема имеет следующий вид [3, 6]:
(—п+1 _ —га)
—-- + аН2) №+1 + (1 - — = щ- (7)
т
Здесь —П = (^П, иХ н, иу н, Ту,) — сеточная вектор-функция решения на п-м временном слое в узле (г, ]), т — шаг по времени, 6 — весовой параметр. Оператор АН2) включает симметричные аппроксимации со вторым порядком первых и вторых пространственных производных на трехточечном шаблоне по каждой пространственной координате. Оператор И,у в правой части равенства (7) состоит из симметричных по каждой пространственной координате аппроксимаций со вторым порядком смешанных производных из уравнений импульсов и слагаемых диссипативной функции из уравнения энергии. Определенная таким образом разностная схема (7) аппроксимирует исходную дифференциальную систему (1) с порядком 0(т+к2) и абсолютно устойчива при 6 > 0.5.
В соответствии с [6] схема (7) после приближенной факторизации реализовывалась методом дробных шагов. По периодической координате х выполнялась циклическая прогонка, алгоритм реализации которой был взят из [7]. В расчетах использовались шаги по времени и пространству к = 0.1 и т = 0.01 соответственно.
2.2. Влияние объемной вязкости на кинематику завихренности
Все расчеты проводились при М = 0.2... 0.5 и И,е = 100. Значение числа Прандтля Рг = 0.74 и отношений 7 = 1.4, а = 0... 2 соответствовали случаю двухатомных газов. При этих значениях параметров диссипативный эффект объемной вязкости рассматривался на модельной задаче в [1, 3].
В сжимаемом нестационарном течении эволюция вихревого возмущения прослеживалась по поведению поля изолиний полной завихренности. Сеточная функция полной завихренности рассчитывалась с помощью аппроксимации со вторым порядком транс-версальной компоненты вектора вихря центральными разностями:
п _ иу,г+1 ,3 ,г— 1 ,3 их,¿,.7+1 их,г,3—1
(■ ■ - — -.
гз 2к 2к '
где их,г,з, иу,г,з — сеточные функции компонент вектора скорости потока, рассчитываемые на основе разностной аппроксимации (7).
В расчетах детально воспроизводилась известная картина динамики крупной вихревой структуры, характерная для возникновения и развития неустойчивости Кельвина— Гельмгольца [5, 8, 9]. В начальный момент времени Ь = 0 изолинии поля завихренности ( для чисел Маха М = 0.2, 0.5 и числа Рейнольдса Ке = 100 представлены на рис. 1. Видно, что вводимое в поток возмущение первоначально сосредоточено в узкой зоне малых скоростей в окрестности у = 0, где возникают замкнутые изолинии, выделяющие пару вихрей равной интенсивности. При этом вне полосы у = ±2 возмущение
исчезающе мало. Из сравнения рис. 1, а и б следует, что в отличие от случая неоднородной несжимаемой жидкости [8, 9], где более коротким волнам соответствуют более интенсивные возмущения, здесь из-за влияния сжимаемости начальные возмущения возрастают с увеличением числа Маха, т. е. с увеличением длины волны.
В последующие моменты времени средняя часть возмущенной области быстро расширяется, выделенные вихри объединяются и образуется близкая к круговой подобласть концентрированной завихренности с замкнутыми изолиниями — так называемое ядро (core). Одновременно по обе стороны от него формируются относительно узкие слои — шнуры (braids), вдоль которых завихренность перетекает с периферии в ядро. В свою очередь, ядро создает в шнурах поддерживающие их растягивающие напряжения. С учетом периодичности в пространстве возникает структура cat's-eye, в которой геометрия изолиний в ядрах топологически эквивалентна особым точкам типа "центр", а на границах раздела — особым точкам типа "седло" [10]. Структура достаточно быстро, при t ~ 2, достигает максимального размера, после чего возмущение затухает. Характерная картина изолиний завихренности, близкая к достигаемому максимуму, показана на рис. 2. Можно заметить, что в этот момент радиус ядра в полтора раза больше, а ширина слоя на границе в полтора раза меньше, чем поперечный размер начальной зоны возмущения на рис. 1. Затухание возмущения продолжается до момента времени t ~ 3.5, после чего структура стабилизируется.
Рис. 1. Картина изолиний поля завихренности О в момент времени £ = 0: а — М = 0.2, Ие = 100; б — М = 0.5, Ие = 100
Рис. 2. Картина изолиний поля завихренности и в момент времени £ = 2.5 для И,е = 100 и М = 0.5: а — а = 0; б — а = 2
Как известно [5], механизм развития неустойчивости Кельвина—Гельмгольца носит невязкий инерционный характер, и влияние динамической вязкости (числа Рейнольдса) на кинематику развития крупной вихревой структуры относительно невелико [8, 9]. В этой связи можно было ожидать, что дополнительный диссипативный эффект объемной вязкости будет проявляться еще слабее.
Степень влияния объемной вязкости на кинематику завихренности оценивалась по зависимостям от времени некоторой условной площади вихревой структуры для различных значений коэффициента а. За границу площади принималась первая незамкнутая изолиния завихренности, входящая в граничную особую точку — "седло". Ей соответствовало значение завихренности ш = -0.15 (см. рис. 2).
Для нахождения временной зависимости условной площади структуры Ба) в расчетной области вводилась дополнительная регулярная сетка с шагами
Дж, = Ду = к = 0.05
так, что центры ее ячеек совпадали с узлами исходной сетки, на которой рассчитывались поля гидродинамических величин. Значения завихренности в дополнительных узлах находились линейной интерполяцией. Искомая площадь определялась следующим образом:
N N
Б(*, а) = ^][](ДжгДу(у),
г=1 з = 1
где Zij имеет вид
ij
uij < -0.15, uij > -0.15.
На рис. 3 приведены графические зависимости величины
п(г,а) = Б (г, а)/Б(г, 0)
для чисел Рейнольдса Ке = 100, Маха М = 0.2, 0.5 и различных значений параметра а. Видно, что расслоение кривых по параметру а, связанное с дополнительным
1.20
1.15
1.10
1.05
32-И ~/Ау/ l-fffff
1-М
1.25
1.20
1.15
1.10
1.05
¡/¿fa
4l//fi 3-fffj ти
1-Ш
Рис. 3. Зависимости от времени относительных площадей n(t,a) для Re = 100 и M = 0.2 (а), M = 0.5 (б): 1 - а = 0; 2 - а = 0.5; 3 - а = 1; 4 - а = 1.5; 5 - а = 2
диссипативным эффектом объемной вязкости, невелико. Для всех значений параметра эволюция структуры носит универсальный характер: рост и достижение максимального значения при £ ~ 2, затем спад и стабилизация при £ > 5 к некоторому значению, несколько превышающему начальное. При сравнении рис. 3, а и б можно заметить, что в противоположность линейной теории [4], в рамках которой наибольшие значения у коротковолновых возмущений (табл. 1), на нелинейной стадии эволюции наибольших значений достигают возмущения с большей длиной волны (с большим числом М). Можно предположить, что к этому эффекту, отмеченному для несжимаемой стратифицированной жидкости в [8], здесь, как и для начальных возмущений, добавляется влияние сжимаемости.
Для количественной оценки влияния объемной вязкости параметра а на расплыва-
ние вихревой структуры вычислялись средние по времени относительные отклонения:
©
(а) = О-1 У а) - 1| ■ 100% а = 0 ... 2.
о
Данные расчетов относительных отклонений (а) для чисел Маха М = 0.2, 0.5, Рей-нольдса Ке = 100 и интервала усреднения по времени О = 6 приведены в табл. 2.
Из данных табл. 2 следует, что рост значения параметра а при фиксированных числах Маха приводит к более интенсивному расплыванию вихревой структуры в пространстве, хотя количественное проявление эффекта, как и ожидалось, незначительно. Более интересно отметить, что при фиксированном значении а расплывание увеличивается примерно пропорционально возрастанию числа Маха. Однако узкий диапазон чисел М в проведенных расчетах не позволяет более полно проследить влияние сжимаемости на кинематику завихренности.
Таблица 2. Относительные отклонения es (а), %, для Ие = 100
а 0.0 0.5 1.0 1.5 2.0
М = 0.2 0.0 0.468 0.953 1.595 2.253
М = 0.5 0.0 0.611 1.239 1.952 2.686
Рис. 4. Зависимости | ш(0, у)| для фиксированных моментов времени Ь при Ие = 100, М = 0.5: а - а = 0; б - а = 2; 1 - Ь = 0; 2 - Ь = 1; 3 - Ь = 2.5; 4 - Ь = 4
Таблица 3. Относительные отклонения еш (а), %, для И,е = 100
а 0.0 0.5 1.0 1.5 2.0
М = 0.2 0.0 0.663 1.546 2.592 3.679
М = 0.5 0.0 0.949 2.216 3.704 5.227
Количественная картина эволюции завихренности во времени представлена на рис. 4, где приведены профили | ш | в поперечном сечении х = 0, проходящем через центр структуры. Известно [11], что в рамках уравнений Навье—Стокса полная циркуляция локализованной вихревой структуры сохраняется во времени. Поэтому о расплывании вихря в пространстве можно судить по уменьшению завихренности в его центре. В данном случае центр структуры совпадает с началом координат, которое, в отличие от разделительной изолинии, ограничивающей площадь Б(г, а), однозначно зафиксировано на расчетной сетке.
Для оценки влияния параметра а на изменения модуля завихренности в центре структуры
| ш(0, г, а)| = ш0(Г, а)
вычислялись относительные отклонения
©
■ 100% ¿г, а = 0 ... 2, о
где интервал усреднения по времени в = 6 можно принять за время "жизни" структуры. Результаты расчетов относительных отклонений еш (а) для чисел Рейнольдса Ке =100 и Маха М = 0.2, 0.5 приведены в табл. 3.
Можно констатировать, что данные табл. 3 с большей определенностью подтверждают выводы о характере влияния параметра а и числа Маха на кинематику вихревой структуры, сделанные на основании табл. 2.
£ш(а) = в
-1
шо(г, а) - Шо (г, 0)
Шо (г, 0)
2.3. Диссипация кинетической энергии возмущений
Наиболее важным с физической точки зрения является влияние объемной вязкости на энергообмен между несущим потоком и возмущениями, который в данной постановке воспроизводится более реально, чем в модельной задаче [1]. Аналогично работам [3, 4] рассматривалась эволюция во времени кинетической энергии возмущения
хо Уо
Е(г) = 11 ¿х I ¿у в «2 + О (8)
-хо -Уо
и абсолютной величины рейнольдсовых напряжений
хо Уо
^ху (г) = J ¿х J ¿у | ви'хи'у | . (9)
- х о - Уо
Величина производства энергии возмущений рассчитывалась на основе полученного в [1] интегрального уравнения энергетического баланса:
¿Е 1
Я(г) = — = 71 + /2 - — (Л + а^), (10)
аг Ке
где слагаемое
хо уо
¿и,
3\ = — J ¿х ! ¿у ди'хи'у-^ (11)
-хо -уо
описывает обмен энергией между возмущением и основным потоком. Интеграл
хо Уо '
* = / ^¡¿ур!^X + дуу) (12)
-хо -уо
интерпретируется как работа при пульсационном расширении (сжатии) газа. Интегралы
хо Уо г 2 2 2 2
ди х\ / +(д^Л +(д^Л +1( д< + диУ
= J ¿х J ¿у
-хо -уо
дх ) \ ду ) \дх ) \ ду ) 3 \ дх ду
(13)
хо Уо 2
*=/ */*(£+а (14)
-хо -уо
определяют вклады диссипативных процессов, связанных со сдвиговой и объемной вяз-костями соответственно. В формулах (11)—(14) величины и'х,и'у,р' — пульсации компонент скорости и давления, хо = п/в, у0 = 10, а волновые числа в для различных чисел Маха М даны в табл. 1. Интегралы вычислялись по квадратурным формулам трапеций с шагом к = 0.1 на использованной в расчетах сетке.
Пульсационные характеристики течения определялись как разности
Ф' (*) = Ф(*) — Ф^),
где компоненты вектор-функции Ф(£) представляют собой мгновенные значения характеристик возмущенного течения, а Ф,(£) — соответствующие характеристики несущего потока. Так как в данном случае, в отличие от модельной задачи [1], невозмущенное течение не является точным стационарным решением системы Навье—Стокса (1), его мгновенные характеристики рассчитывались параллельно с расчетом возмущенного потока.
Для контроля вычислений значений производства энергии 0(£) из энергетического уравнения (10) значения 0(£) рассчитывались непосредственно с использованием (8) и симметричной конечно-разностной аппроксимации временной производной:
о® = Е " + т >-Е — т).
Как показали расчеты, альтернативные значения 0(£) различались не более, чем на 1.5%.
Временные зависимости кинетической энергии возмущений Е(¿, а) представлены на рис. 5. Чем больше значение параметра а, тем меньше кинетическая энергия структуры на всем временном интервале. Эволюция Е(¿, а) соответствует кинематике развития структуры, намеченной в комментарии к рис. 2. С ростом структуры возрастает пуль-сационная энергия, которая после достижения максимума спадает и стабилизируется к некоторому значению на временах в > 6. Сдвиг максимумов на графиках рис. 2 и 5 связан с условностью вычисления величины площади структуры Б(¿, а), в то время как
2
а б
Рис. 5. Временные зависимости энергии возмущений для Ие = 100 и М = 0.2 (а),
М = 0.5 (б): 1 - а = 0; 2 - а = 0.5; 3 - а = 1; 4 - а = 1.5; 5 - а = 2
энергия Е(г, а) вычисляется интегрированием по всей расчетной области. Временные зависимости модуля рейнольдсовых напряжений аху(г, а), определенных в (9), практически повторяют приведенные графики Е(г, а).
Для количественного сравнения вклада объемной вязкости в диссипацию кинетической энергии возмущений вычислялись средние по времени относительные отклонения
©
Ее (а) = 0-1у
о
Осреднение проводилось по условному времени "жизни" структуры 0 = 6. Результаты расчетов относительных отклонений Ее (а) представлены в табл. 4. Из табл. 4 следует, что рост параметра а при фиксированном числе Маха приводит к большей диссипации кинетической энергии возмущения. При числе Маха М = 0.5 среднее по времени относительное уменьшение энергии для а = 2 достигает Ее = 13.567%. Это хорошо согласуется с аналогичными данными, полученными в [1, 3] для модельной задачи, хотя в ней воспроизводилось только затухание структуры.
Расчетные зависимости производства пульсационной энергии Д (г), а также отдельных вкладов .]1, .а, .]4 в общий баланс процесса, полученные из энергетического уравнения (10), даны на рис. 6-9. Анализ графиков на этих рисунках позволяет детализировать процесс эволюции кинетической энергии структуры и ее диссипации. На временном интервале, где происходит рост структуры и кинетической энергии возмущения, ее производство сначала положительно и возрастает, достигая некоторого максимума, а затем начинает убывать и в конце концов становится отрицательным, что соответствует убыванию энергии структуры. Имеет место точное соответствие характерных
Е(г, а) - Е(г, 0)
Е(г, 0)
100% ¿г, а = 0 ... 2.
Таблица 4. Относительные отклонения ее (а), %, для Ие = 100
а 0.0 0.5 1.0 1.5 2.0
М = 0.2 0.0 1.347 3.148 5.287 7.511
М = 0.5 0.0 2.478 5.772 9.633 13.567
точек на графиках рис. 5 и 6. В частности, точка перехода кривых через нуль соответствует максимуму на кривых Е(¿), а точки максимума и минимума отвечают точкам перегиба на графиках Е(¿).
Сравнение графиков рис. 6 и 7 показывает, что производство энергии возмущений почти полностью определяется вкладом интеграла J1 (¿), характеризующим процесс обмена энергией между возмущением и несущим потоком в уравнении (10). В подынтегральном выражении Jl(t) производная от скорости несущего потока бПа/бу > 0 на всем расчетном интервале времени. Вместе с тем, на первом этапе эволюции структуры, когда J1 (t) > 0, рейнольдсовы напряжения дп'хи'у также, по крайней мере, неотрицательны в большей части расчетной области. При этом происходит подкачка энергии в вихревую структуру из внешнего течения, что и определяет рост ее кинетической энергии. Такое направление потока энергии в феноменологической теории турбулентности принято связывать с явлением "отрицательной" вязкости. Переход Jl(t) в отрицательную
Рис. 6. Временные зависимости производства пульсационной энергии для И,е = 100 и
М = 0.2 (а), М = 0.5 (б): 1 — а = 0; 2 — а = 0.5; 3 — а = 1; 4 — а = 1.5; 5 — а = 2.
а б
Рис. 7. Временные зависимости интеграла для И,е = 100 и М = 0.2 (а), М = 0.5 (б): 1 — а = 0; 2 — а = 0.5; 3 — а = 1; 4 — а = 1.5; 5 — а = 2
область определяется сменой знака рейнольдсовых напряжений, которые становятся отрицательными в большей части области и создают обычный для них диссипативный эффект, способствуя затуханию возмущения.
Сопоставляя графики рис. 6-9, можно заключить, что отрицательный вклад дисси-пативных слагаемых Дз(£) и Д4(£) оказывает более ощутимое влияние на производство кинетической энергии возмущения на временах, близких к положительным максимумам и Д\(Ь), что видно из сравнения графиков рис. 6 и 7. Вместе с тем можно констатировать, что нули этих графиков имеют незначительный сдвиг по времени, а отрицательное производство в области О(Ь) < 0 связано в основном с поведением Д\(1) < 0. В частности, совпадает характер расслоения обеих зависимостей по параметру а, в которых большим значениям объемной вязкости соответствуют меньшие скорости затухания О(Ь) (менее отрицательные значения <7*1 (¿)). Такое поведение противополож-
Рис. 8. Временные зависимости интеграла Jз(t) для Ие = 100 и М = 0.2 (а), М = 0.5 (б): 1 - а = 0; 2 - а = 0.5; 3 - а = 1; 4 - а = 1.5; 5 - а = 2
12 34 5* 12 345/
а б
Рис. 9. Временные зависимости интеграла J4(t) для Ие = 100 и М = 0.2 (а), М = 0.5 (б): 1 - а = 0; 2 - а = 0.5; 3 - а = 1; 4 - а = 1.5; 5 - а = 2
Таблица 5. Эволюция во времени интеграла ^(а) для И,е = 100, М = 0.5
а 0.0 0.5 1.0 1.5 2.0
г = 1 0.036219 0.036590 0.036144 0.035265 0.034306
г = 2 -0.039891 -0.037281 -0.034514 -0.031827 -0.029524
г = з 0.000170 0.000300 0.000270 0.000213 0.000197
г = 4 0.021402 0.020719 0.020008 0.019332 0.018742
г=5 -0.005725 -0.005292 -0.005101 -0.005091 -0.005002
но результату, полученному в [1], где скорость затухания возрастала пропорционально увеличению а. Это означает, что при Д(г) < 0 затухание кинетической энергии возмущения не определяется отрицательными вкладами Л3(г), Л4(Г). Следовательно, расчетное число Рейнольдса Ке = 100 значительно превышает критическое значение числа Рейнольдса, которое находится из решения вариационной задачи энергетической теории устойчивости [11], получаемой из уравнения (10) при Д(г) = 0.
Сравнение друг с другом значений интегралов Л3(Г) и Л4(Г), приведенных на рис. 8 и 9, показывает, что при М = 0.2 диссипативный вклад сдвиговой вязкости на порядок превышает соответствующий вклад объемной вязкости. Однако уже для М = 0.5 оба вклада становятся сравнимыми Л3(Г) ~ Л4(Г), особенно с учетом а > 1. Это позволяет предположить, что при более сильном возбуждении, затрагивающем колебательные степени свободы, диссипативный эффект, связанный с релаксационным процессом, может стать определяющим фактором подавления возмущений, как это было получено в [2] на модельной задаче.
Выборочные значения интеграла Л2(Г,а), описывающего процессы пульсационного расширения и сжатия в уравнении (10), в различные моменты времени Г для М = 0.5 и Ке = 100 даны в табл. 5.
Сравнивая данные табл. 5 с характерной величиной остальных интегралов Л1, Л3, Л4 энергетического баланса, можно констатировать, что они разнятся более чем на два порядка и вклад Л2(Г) в первом приближении пренебрежимо мал. Рассмотрение подынтегральных выражений в Л2, Л4 показывает, что это может быть объяснено малостью пульсаций давления р', которые через уравнение состояния (см. систему уравнений (1))
связаны с пульсациями температуры Т' и плотности в' соотношением
р' = (вТ' + в'Тв).
В работе [12] уравнение энергетического баланса было записано относительно производства полной энергии возмущений, которая выражалась интегралом
™ = / { ^ + [(£ ) +<7-1- (С)]} а".
П
Здесь дополнительные слагаемые по сравнению с интегралом кинетической энергии Е(Г) в (10) получены из преобразования интеграла Л2(Г). Сравнение величины отдельных вкладов показывает, что для определенного класса сжимаемых течений можно в энергетической теории устойчивости пренебречь интегралом Л2 (Г) в уравнении вида (10) и ограничиться рассмотрением производства только кинетической энергии возмущений.
Выводы
Выполнено численное моделирование нелинейного развития дозвуковых вихревых возмущений в эволюционирующем во времени сдвиговом слое. Диапазон параметров течения соответствовал реальным значениям для двухатомных газов. Полученные результаты позволили количественно оценить влияние объемной вязкости на формирование и эволюцию развивающейся естественным образом вихревой структуры и сделать следующие выводы.
1. Влияние объемной вязкости на кинематику вихревой структуры с преимущественно невязким механизмом эволюции относительно мало, но ее дополнительный диссипа-тивный эффект проявляется вполне определенно.
2. При возрастании объемной вязкости в реальных для двухатомных газов пределах относительное увеличение диссипации кинетической энергии структуры, осредненное по времени ее "жизни", достигает 13 %, что хорошо согласуется с аналогичными результатами предыдущих исследований на модельной задаче [1, 3].
3. В дозвуковом диапазоне чисел Маха при M > 0.5 вклады объемной вязкости (сжимаемости) и сдвиговой вязкости в темп диссипации кинетической энергии возмущений становятся соизмеримыми по порядку величины.
4. Для возмущенных плоскопараллельных дозвуковых течений эволюция энергии возмущений с достаточной точностью описывается функционалом кинетической энергии пульсаций.
Список литературы
[1] Григорьев Ю.Н., Ершов И.В. Подавление вихревых возмущений релаксационным процессом в течениях возбужденного молекулярного газа // ПМТФ. 2003. Т. 44, № 4. С. 22-34.
[2] Григорьев Ю.Н., Ершов И.В., Ершова Е.Е. Влияние колебательной релаксации на пульсационную активность в течениях возбужденного двухатомного газа // ПМТФ. 2004. Т. 45, № 3. С. 15-23.
[3] Григорьев Ю.Н., Ершов И.В., Зырянов К.И., Синяя А.В. Численное моделирование эффекта объемной вязкости на последовательности вложенных сеток // Вычисл. технологии. 2006. Т. 11, № 3. С. 36-49.
[4] Blumen W. Shear layer instability of an inviscid compressible fluid //J. Fluid Mech. 1970. Vol. 40, part 4. P. 215-239.
[5] Бетчов Р., Криминале В. Вопросы гидродинамической устойчивости / Под ред. О.Ф. Васильева и В.В.Пухначева. М.: Мир, 1971.
[6] Ковеня В.М., Яненко Н.Н. Метод расщепления в задачах газовой динамики. Новосибирск: Наука. Сиб. отд-ние, 1981.
[7] Квасов Б.И. Интерполяция кубическими и бикубическими сплайнами: учеб. пособие. Новосибирск: НГУ, 2004.
[8] Patnaik P.C., Sherman F.S., Corcos G.M. A numerical simulation of Kelvin—Helmholtz waves of finite amplitude //J. Fluid Mech. 1976. Vol. 73, part 2. P. 215-239.
[9] Corcos G.M., Sherman F.S., The mixing layer: deterministic models of a turbulent flow. Part I. Introduction and two-dimensional flow //J. Fluid Mech. 1984. Vol. 139. P. 29-65.
[10] Степанов В.В. Курс дифференциальных уравнений. М.: ФМ, 1959.
[11] Ламб Г. Гидродинамика. М.; Л.: ОГИЗ—Гостехихдат, 1947.
[12] Grigoryey Yu.N., Ershoy I.V. Influence of bulk viscosity on the Kelvin—Helmholts instability // Proc. XIII Intern. Conference on the Methods of Aerophys. Res. 2007. Part III. P. 123-128.
Поступила в редакцию 22 февраля 2008 г., в переработанном виде — 3 марта 2008 г.
