CC BY
25
Вычислительные технологии
Том 15, № 6, 2010
Численное моделирование инерционной неустойчивости в колебательно неравновесном
двухатомном газе*
Ю.Н. Григорьев Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: grigor@ict.nsc.ru
И. В. Ершов, К. И. Зырянов Новосибирский государственный архитектурно-строительный
университет, Россия e-mail: grigor@ict.nsc.ru
В рамках уравнений двухтемпературной гидродинамики исследовано влияние колебательной релаксации на подавление неустойчивости Кельвина—Гельмгольца в развивающемся во времени сдвиговом слое колебательно неравновесного двухатомного газа.
Ключевые слова: неустойчивость Кельвина—Гельмгольца, колебательная релаксация, кинетическая энергия возмущений, диссипация.
Введение
В работе [1] рассматривалось влияние умеренного возбуждения внутренних степеней свободы молекул на развитие инерционной неустойчивости Кельвина—Гельмгольца в свободном слое сдвига в молекулярном газе. Расчеты были выполнены на основе полных уравнений Навье—Стокса сжимаемого теплопроводного газа, в которых с помощью коэффициента объемной вязкости учитывается только возбуждение вращательных мод, а колебательные степени свободы предполагаются "замороженными". Несмотря на инерционный характер возбуждения волн Кельвина—Гельмгольца результаты расчетов показали, что с возрастанием объемной вязкости усиливается диссипация кинетической энергии пульсаций, а скорость ее производства уменьшается. Вместе с тем было высказано предположение, что при дополнительном возбуждении колебательных уровней молекул данные эффекты должны возрастать. В этом случае релаксационный процесс может быть использован для управления течением, например, для затягивания ламинарно-турбулентного перехода, так как колебательная неравновесность легко создается искусственным путем.
В данной работе нелинейное развитие неустойчивости Кельвина—Гельмгольца рассматривается в рамках уравнений двухтемпературной газовой динамики. В них предполагается, что поступательные и вращательные степени свободы образуют квазиравновесный термостат, характеризуемый статической температурой потока, а релаксация колебательных степеней к равновесию описывается уравнением Ландау—Теллера для
*Работа выполнена при финансовой поддержке РФФИ (грант № 08-01-00116).
колебательной температуры газа [2, 3]. В качестве начальных возмущений использовались невязкие вихревые моды с максимальным инкрементом нарастания. Последние предварительно рассчитывались на основе линеаризованной системы уравнений невязкой двухтемпературной газовой динамики.
1. Постановка задачи
1.1. Основные уравнения и начально-краевые условия
В координатной плоскости (x1;x2) рассматривается плоскопараллельное сдвиговое течение двухатомного газа, в котором стационарный (несущий) поток равномерной плотности р0 направлен вдоль оси абсцисс x1 и аналогично [1] имеет профиль скорости с точкой перегиба при x2 = 0
Us(x2) = Uo • th(x2/^o)-
Считается, что в невозмущенном течении все внутренние степени свободы молекул находятся в равновесии и описываются единой статической температурой потока T0 = const, В отсутствии внешней накачки энергии в колебательные степени свободы это единственно возможное стационарное решение энергетических уравнений используемой системы. Характерная длина для обезразмеривания определялась параметром S0, описывающим максимальный наклон профиля скорости
So = Uo
'dUa
dx2 Х2=0_
-1
В качестве других масштабирующих величин были выбраны асимптотическое значение скорости [/о, постоянные плотноеть р0 и температура Т0 и образованные из них характерные время т0 = 50/и0 и давление р0 = р0ио2. В определенных на их основе безразмерных переменных система уравнений двухтемпературной гидродинамики имеет вид
dp Opui dt dxi
0,
( дщ dt
dui
р\~ж+и^
dp 1 d2ui dxi Re dx2
1
1\ d2 u,
^ Re V 1 ^ 3 ) dxi dxj '
(dT
dT
Y
Re Pr dx2
d2T РЪгь{ТугЬ~Т) |
y(y - 1)m2
2Re
dui du.
+ 2 ai -
2
Tvib (1 — Yvib) 2
duu
3 / \ dx.
I dTvib
p[-dr
u,
dTvib\ Yo.2 (1 — Yvib) d2Tvib 0 (T —
dxj
Yvib Re Pr dx2
Tvib
YM2p = pT, i,j = 1, 2.
(1)
Здесь и далее в записи уравнений в координатной форме по повторяющимся индексам подразумевается суммирование.
Параметры, входящие в уравнения системы (1), определяются следующим образом: а\ = пь/п ~ отношение объемной вязкости п> к сдвиговой вяз кости п; а2 = + Ког), где Кы, Ки> — коэффициенты теплопроводности, описывающие
соответственно упругие энергообмены между поступательными степенями свободы и неупругие энергообмены вращательных и колебательных степеней свободы с поступательными; Tvib — характерное время релаксации колебательных мод молекул к равновесному состоянию; коэффициент y = cp/cv - показатель адиабаты, где Cv = "cv.tr+Cv, rot, cp = cv + R — соответственно удельные теплоемкости при постоянных объеме и давлении, связанные с поступательным и вращательным движением молекул, R — газовая постоянная; коэффициент Yvib = cv,vib/(cv + cv,vib) характеризует степень неравновесности колебательных мод, cv, vib — парциальная удельная теплоемкость при постоянном объеме, характеризующая запас энергии в системе возбужденных колебательных мод; параметры Ее = Uoóopo/v и М0 = Uo/y/jRTo — числа Рейнольдса и Маха несущего потока; Pr = ncp/A0 — число Прандтля, Предполагается, что в (1) все коэффициенты переноса и удельные теплоемкости постоянны. Из соотношений Эйкена [2]
л _ 5// Cv> tr л _ 6 7] Cv> rot , _ 6 7] Cv> vib
Л-tr ~ j ^rot Z j ^vib Z
2 5 5
следует удобное для использования в расчетах выражение
Avib 20 Tvib
Atr + Arot 33 (1 — Yvib) В безразмерных переменных стационарный несущий поток задается соотношениями
Us{x2) = thx2, Ts = Tvib>s = ps = 1, ps = (2)
Исходная краевая задача ставится в бесконечной полосе, центр которой совпадает с началом координат: x1 G [—xi,0, xi,0], x2 G (—о, оо). Ширина полосы по координате x1 выбиралась равной длине волны возмущения A = 2п/в При этом x1,0 = п/в здесь в — волновое число, В расчетах асимптотические условия при x2 — ±о переносились на x2 = ±x2,0, где ордината x2,0 определялась из условия | thx2,0 — l| < 10-12, В итоге было принято значение x2,0 = 20,
В момент времени t = 0 на основной поток накладывалось двумерное возмущение с длиной волны A и волновым вектором к = (в, 0), Начальные условия для поля скорости и термодинамических величин, включающие возмущения, определялись в виде
u1(0, x1, x2) = thx2 + u1(0, x1, x2), u2(0, x1, x2) = u2(0, x1, x2),
T(0, X1, X2) = 1 + T'(0, X1, X2), Tvib(0, X1, X2) = 1 + T'ib(0, X1, X2),
p(0, xi, x2) = 1 + p'(0, xi, x2), p(0, xi, x2) =+p'(0, xi, x2). (3)
В качестве вводимых в основной поток начальных возмущений компонент вектора скорости «1, u2 и термодинамических величин р TT'ib, p' использовались собственные линейные невязкие колебания с наибольшими инкрементами нарастания.
При t > 0 на границах x1 = ±x1,0, x2 G [—x2,0,x2,0] ставились периодические условия
«1(t, Х1, 0, X2) = U1(t, —X1, 0, X2), «2(t, X1, 0, X2) = «2(t, —X1, 0, X2),
T (t, X1, 0, X2) = T (t, —X1, 0, X2), Tvib(t, X1, 0, X2) = Tvib(t, —X1, 0, X2),
p(t, Xi, 0, Х2) = p(t, -Xi, 0, X2), p(t, Xi, 0, X2) = p(t, -Xi, 0, X2), a при x2 = ±X2, 0, Xi G [—Xi, 0, Xi, 0] принимались условия невозмущенного потока
Ui(t, Xi, X2, 0) = -Ui(t, Xi, -X2, 0) = 1, U2(t, Xi, X2, 0) = U2(t, Xi, -X2, 0) = 0,
T(t, Xi, X2, 0) = T(t, Xi, -X2, 0) = 1J Tvib(t, Xi, X2, 0) = Tvib(t, Xi, -X2, 0) = 1,
1
p(t, Xi, X2, 0) = p(t, Xi, -X2, 0) = 1, p(t, Xi, X2, 0) = p(t, Xi, -X2, 0)
Y M2'
1.2. Расчет начальных возмущений
Параметры собственных линейных невязких колебаний с наибольшими коэффициентами нарастания определяются как решения системы уравнений двухтемпературной газовой динамики, линеаризованной относительно стационарного основного течения (2):
+ u ■— + — = о ди* I и ди'1 I др' = о
dt Я,г ÔXi ÔXi ' dt s,j dXj ÔXi
— + u. dT'\h i)d< I =0
dt s'3 dxj dxi Tvib(l — Yvib)
dT' dT' (T ' — T' )
vib t T T vib \± ± vib ) r\ ]\ /Т ^ ' / , r-p' • • 1 r, ( i\
— ' ' rmb =0' ^MP=P+T> ^ = 2- (4)
В отличие от случаев идеального газа [4] и идеальной жидкости [5] соответствующие решения системы (4) ранее не рассматривались, С использованием стандартной схемы [4, 5] искомые решения представляются в виде плоских волн, распространяющихся в направлении основного потока:
{u/i,u/2,p',T',T'mb,p'} = {u(x2),v(x2),p(x2),9(x2),9Vib(x2),p(x2)}exp [0(Xi - et)], (5)
где u, v, p, d, dvib, p — амплитудные функции возмущений, i — мнимая единица, c = cr + ici — комплексная фазовая скорость волны.
Подстановка (5) в систему (4) приводит к спектральной задаче, которая при cr = 0 была преобразована к скалярному уравнению для амплитуды пульсаций давления p:
d 2p 2 sech2x2 dp
T~2 ~ 7TÛ-1--P2 1 ~ m2 (th ж2 - ici)
dx2 (th x2 — ici) dx2 L 4 '
p = 0,
pL = -œ= p\x2=+œ= 0, p(x2) = pr (x2) + ipi (x2) ,
I 1 + Yv + ifîTvib (th X2 - ici) Yvib fn,
™ = M - x -—-—, 7„ = --, 6)
1 + (Yv/y) + ipTvib (th X2 - ici) 1 - Yvib
где pr, pi — действительная и мнимая части комплексной собственной функции давления p, а собственными значениями задачи (6) являются q. Параметрами задачи служат волновые числа в числа M axa M коэффициент Yvib, характеризующий степень возмущения колебательной моды, и время колебательной релаксации rvib.
При данном профиле скорости Us\x2^±^= ±1 и US\x 0 и асимптотика урав-
нения (6) при х2 — ±œ имеет вид
p''- в2 [1 - m2±)M0(1 T Q)2] p = 0.
Отсюда асимптотическое решение определяется как
Р ~ ехр (тРх2^1-т2±)М2(1Тсг)2) . (7)
2 _ 1 + Ъ ± iprvib{ 1 =F Ci)
Д6СЬ Ш(±) " 1 + (ъИ)±гРттЬ(1тсгУ Для численных расчетов собственных значений неустойчивых невязких мод использовалась методология работ [4, 5]. Вводились новые переменная z = th x2 и функция
„ ч d ln p
П(ж2) = -Г^ = Пг + Шг,
dx 2
в которых уравнение (6) принимает вид
£Ш _ /32[1 — m2(z — icj)2] — П2 | 2П dz (1 — z2) (z — ici) '
nL=±i= ТР-^1~т2{±)(1тгсг)2, -1<z<1. (8)
Граничные условия в задаче (8) находятся из асимптотики (7), Для расчетов методом "стрельбы" необходимы также значения производных функции П на концах интервала интегрирования. Они определяются непосредственно из уравнения (8) при z — ±1, Возникающая при этом в первом слагаемом правой части (8) неопределенность типа [0/0] раскрывается по правилу Лопиталя, В результате имеют место выражения
dU 132(1тгсг)2- 2 m2±) Т K(±) (1 Т ici) -4/3- ^1-т2±)(1тгсг)2
dz z=±1 2 (1 т ici) ■ 1 + /3-^ т2±) (1 Т ici)2
ЪгЪ ФТъгЬ^ - 1)М2 Л(±) = —7--Г ' -о- (У)
7(1 - Ъгъ) [1 + (7„й/(7(1 - 7^гъ)) ± гвттъ (1 т гсг)]
Для расчета спектральных характеристик неустойчивых мод из уравнения (7) выделялись уравнения для действительной Пг и мнимой П» частей функции П, Полученная таким образом система при фиксированных наборах параметров интегрировалась с данными Коши (8), (9) с помощью процедуры Рунге—Кутты четвертого порядка на интервалах г € [-1, 0] и £ € [0, 1] с шагом Дz = 10-3, Волновые числа изменялись в пределах в € [0, 1] с шагом Дв = 0.1.
Точкой "прицеливания" служило £ = 0, Значения с подбирались таким образом, чтобы вычисленные "слева" и "справа" в точке £ = 0 значения функций Пг, П» совпадали с точностью до 10-8, Соответствующее такому совпадению значение с принималось в качестве собственного при заданном наборе параметров, В табл. 1 приведены примеры найденных собственных значений фазовых скоростей с», соответствующих волновых чисел в и инкрементов роста вс» наиболее неустойчивых невязких мод для несущего потока с числами Маха М = 0.2, М = 0.5, временами колебательной релакеации т^ъ = 0.1, т^ъ = 1 и значениями параметра 7^ъ = 0 — 0.4,
После нахождения максимальных инкрементов нарастания рассчитывались соответствующие им собственные функции спектральной задачи (6), Уравнение (6) для
Таблица 1. Спектральные характеристики и инкременты роста наиболее растущих невязких возмущений для т^ъ = 0.1 и т^ъ = 1
М = 0.2 М = 0.5
ЪгЬ Р Сг Рсг Р Сг Рсг
ТтЬ = 0.1
0.0 0.4260 0.425523 0.181306 0.3970 0.355576 0.141332
0.1 0.4386 0.412995 0.181140 0.3950 0.353816 0.139757
0.2 0.4384 0.412453 0.180819 0.3930 0.351167 0.138009
0.4 0.4379 0.410737 0.179862 0.3849 0.345069 0.132817
ТугЬ 1
0.0 0.4260 0.425523 0.181306 0.3970 0.355576 0.141332
0.1 0.4387 0.412997 0.181182 0.3951 0.354300 0.139984
0.2 0.4385 0.412577 0.180915 0.3932 0.352276 0.138515
0.4 0.4377 0.411503 0.180115 0.3890 0.344847 0.134145
Р = Рг + Рг разделялось на действительную и мнимую части, которые приводились к нормальной системе уравнений первого порядка:
/' = В/,
где f = (рЗг,Рг,Рг,Рг),
(
в
0 0
V
в2(1 - А)
-в 2Аг
0 0
в 2Аг в2(1 - Аг)
1 0
2Ш х28ес112х2
(Ш2 ж2 + с2) 2сг 8ес112х2
0 1
2сг 8ес112х2
(Ш2 ж2 + с2) 2Ш х28ес112х2
(Ш2Ж2 + С2) (Ш2 Х2 + с2) /
Аг = т^Ш2 х2 — с2) + 2сгт2Ш х2, = т2(Ш2 х2 — с2) — 2сгт2^Ш х2,
тг(х2) и тг(х2) — соответственно вещественная и мнимая части функции ш(х2), определенной в (6), Полученная система интегрировалась с использованием процедуры Рунге—Кутты четвертого порядка па интервалах х2 € [0, —20] и х2 € [0, 20] с шагом Дх2 = 10-3.
х2 = 0
цнп определяются с точностью до постоянного множителя, то нормировка функций р была выбрана из условия, чтобы при — 0 и М = 0 они совпадали с собственными функциями возмущений р для идеального газа [4], В результате начальные условия при х2 = 0
Рг(0) = в2, Рг(0) = 0, рГ(0) = 0, р г(0) = в2Пг(0).
Другие амплитудные функции определялись через найденные функции Р и ¿Р 'с помощью следующих из (4), (5) соотношений:
1 Л _ Л 8всЬ2Х2 Л и = — — грр--——р
В
В
у = -ЪР
р
рм2 (1 + 7у + ттЪР)
1 + (т^/т) + ТтЪ В
1
Рис. 1. Картина изолиний поля завихренности ш в момент времени t = 0 для M = 0.5 и Tvib = 1: а — Yvib = 0\ б — Yvib = 0.4
е = рМ2 (7-1) (1 + ттЬР) е = (7-1)рМ2
1 + (т«/т) + TvibD ' тЪ 1 + {ъ h) + TvibD ' = 7отЬ . I) = г/3 (thx2 - ici).
1 - Yvib
Полученные решения использовались в дальнейших расчетах в начальных условиях (3). Вводимые возмущения определялись как вещественная часть рассчитанных решений (8):
g'(xi, Х2, t) = Re [go(x2) ■ exp(i^xi)] ■ exp(ficrf).
Здесь g' = (u[, u'2, Q T Tvib, p ') — вектор возмущений, g0(x2) = (u, v, g, в, 9vib, p) — вектор соответствующих амплитуд.
На рис. 1 представлены примеры изолиний поля завихренности ш = u2,x — u1,y в начальный момент времени t = 0 для числа Маха M = 0.5, времени колебательной релаксации Tvib = 1 и двух значений степени неравновееноети колебательной энергии 7vib = 0, Yvib = 0.4.
Из сравнения рисунков 1, а, б следует, что возрастание степени неравновесности колебательной энергии молекул газа 7vib приводит к уменьшению амплитуд начальных возмущений. Расчеты также показали, что варьирование значений параметра Tvib от нуля до единицы не оказывает влияния на топологию изолиний и величины амплитуд начальных возмущений.
2. Эволюция возмущений в колебательно неравновесном двухатомном газе
2.1. Разностная схема
Расчеты эволюции возмущений велись на равномерной сетке с шагом к по обеим пространственным переменным. Система (1) аппроксимировалась весовой конечно-разностной схемой с расщеплением по физическим процессам и пространственным переменным, структура которой была заимствована из [6] и уже использовалась в работах [1, 7]. В операторной форме схема имеет вид [6, 7]
(а п+1 — а п)
у-^ + и + =С£, (10)
т
где яп = (рп, иПу, и^уТД^) — сеточная вектор-функция решения па и-м временном слое в узле (г, ]), т — шаг то времени, 6 — весовой параметр. Оператор Ь
включает симметричные аппроксимации со вторым порядком первых и вторых пространственных производных на трехточечном шаблоне по каждой пространственной координате. Оператор G^ составлен из симметричных по каждой пространственной координате аппроксимаций со вторым порядком смешанных производных из уравнений импульсов, первых производных, входящих в диссипативную функцию из уравнения энергии, и нсточннковых слагаемых из релаксационного уравнения и уравнения энергии, В квазиравновесном случае, когда система уравнений релаксационной гидродинамики (1) переходит в систему полных уравнений Навье—Стокса, схема (10) совпадает со схемой, предложенной в [6]. Проведенный анализ (10) показал, что появление релаксационного уравнения и источпикового слагаемого в уравнении энергии не меняет вычислительных свойств схемы по сравнению со схемой [6] для системы Навье—Стокса, Таким образом, разностная схема (10) на регулярной сетке с шагом h аппроксимирует исходную дифференциальную систему (1) с порядком 0(т + h2) и абсолютно устойчива при ö > 0.5,
В соответствии с [6] схема (10) после приближенной факторизации реализовывалась методом дробных шагов. По периодической координате x выполнялась циклическая прогонка, алгоритм реализации которой был взят из [8]. В расчетах использовались шаги по времени и пространству h = 0.1 и т = 0.01,
2.2. Влияние колебательной релаксации на кинематику завихренности
Расчеты проводились при числах Рейпольдса Re = 100, Маха М = 0.2, М = 0.5 и Прандтля Рг = 0.74. Показатель адиабаты был взят 7 = 1.4, Значения степени неравновесности колебательной энергии 7vib, времени колебательной релакеации Tvib и параметра а = п/Пь принимались равными Yvib = 0, 0.1, 0.2, 0.4, Tvib = 0, 0.1, 1 и а = 0, 1, 2, При этих значениях параметров, соответствующих случаю двухатомных газов, днссп-пативный эффект, обусловленный колебательной релаксацией, рассматривался ранее на модельной задаче в [9].
Эволюция вихревого возмущения прослеживалась по поведению поля изолиний полной завихренности. Сеточная функция полной завихренности рассчитывалась с помощью аппроксимации со вторым порядком траневереальной компоненты вектора вихря центральными разностями
<7_ /1 /1 _ /1
п _ 2, г-\-1 j 2, i—lj _ ul,ij+1 ul,ij-l
2h 2h '
где un ij, пП ij — сеточные функции компонент вектора скор ости потока на n-м временном слое.
Расчеты позволили детально воспроизвести известную картину нелинейной динамики крупной вихревой структуры "cat's—eve" в процессе возникновения и развития неустойчивости Кельвина—Гельмгольца [1, 10]. Структура достаточно быстро при t ~ 2.1 достигает максимального размера, затем возмущение затухает до момента времени t & 3.5, после чего структура стабилизируется. Пример изолиний завихренности в структуре вблизи к максимальному размеру при t = 2.5 приведен на рис, 2, Можно заметить, что в этот момент времени радиус ядра в полтора раза больше, а ширина слоя на границе в полтора раза меньше поперечного размера начальной зоны возмущения на рис, 1,
Степень влияния неравновесности колебательной энергии на кинематику завихренности оценивалась по зависимости от времени некоторой условной площади вихревой
3
2
1
0
-1
-2
-3
■4
-6 -4 -2 0 2 4 6 Х1
Рис. 2. Картина изолиний поля завихренности ш в момент времени £ = 2.5 для И,е = 100, М = 0.5, т^ь = 1, а = 0 7™ь = 0.4
структуры для различных значений параметра YmЬ■ За внешнюю границу структуры принималась первая незамкнутая изолиния завихренности, входящая в граничную точку — "седло", которой соответствовало ш = -0.15 (см. рис. 2).
Для нахождения временной зависимости условной площади структуры Б (£, 7^гь) в расчетной области вводилась дополнительная регулярная сетка с шагами Дя1,г = Дхг,^ = к = 0.05 так, что центры ее ячеек совпадали с узлами исходной сетки, на которой рассчитывались поля гидродинамических величин. Значения завихренности в дополнительных узлах находились линейной интерполяцией. Искомая площадь определялась следующим образом:
На рис. 3 приведены примеры графических зависимостей величины п(£, 7™ь) = Б(£, 7тЬ)/Б(£, 0) для различных режимов. Видно, что на участке инерционного нарастания структуры вплоть до £ ~ 2 расслоение кривых по параметру 7^гЬ практически отсутствует. В процессе диссипативпого затухания расслоение также невелико. Для количественной оценки влияния параметра 7^гЬ на раеплывание вихревой структуры вычислялись средние по времени относительные отклонения
где интервал усреднения по времени в = 6 приближенно соответствует времени "жизни" структуры. Данные расчетов относительных отклонений eS(yvib) для ряда режимов приведены в табл. 2. Из таблицы следует, что возрастание параметра Yvib при фиксированных ai, Tvib и числах M axa M приводит к более интенсивному раеплыванию вихревой структуры в пространстве, хотя количественное проявление эффекта незначительно. Расчеты показали, что варьирование параметра Tvib слабо влияет на поведение зависимостей условных площадей S(t, Tvib), Из графиков па рис. 3 и данных табл. 2
©
о
Рис. 3. Зависимости от времени относительных площадей ц(Ь,ЪгЪ) для И,е = 100, М = 0.5, т-югЬ = 1 (шриховая линия — а1 = 0 сплошная — а.1 = 2; 1 — ^тЬ = 0, 2 — ^тЪ = 0.2, 3 — ЪгЬ = 0.4
Таблица 2. Относительные отклонения es (ЪгЬ), % для И,е = 100, т^ь = 1
«1=0 (хх = 2
м ЪгЬ = 0.1 ЪгЬ = 0.2 ЪгЬ = 0.4 ЪгЬ = 0.1 ЪгЬ = 0.2 ЪгЬ = 0.4
0.2 0.046 0.144 0.788 0.056 0.2190 0.841
0.5 0.158 0.326 1.136 0.195 0.372 1.186
следует, что в рассмотренном диапазоне изменения параметров «1 и 7^Ь объемная вязкость оказывает более сильное влияние на раеплывание завихренности, чем релаксация возбужденных колебательных уровней, хотя при более глубоком возбуждении, которое еще можно рассматривать на основе системы уравнений (1), влияние релаксации может возрасти.
2.3. Диссипация кинетической энергии возмущений
Наибольший интерес с физической точки зрения представляет исследование влияния колебательной релаксации на диссипацию кинетической энергии возмущений. В данной постановке задачи этот процесс воспроизводится более реалистично, чем в модельной задаче [9]. Аналогично работам [1, 9] рассматривалась эволюция во времени кинетической энергии возмущения
Х1, 0 Х2, 0
= ^ У (1х\ J (1x2 р ■ + и®2) (11)
—Х1, 0 — Х2, 0
и абсолютной величины рейнольдеовых напряжений
Х1, 0 Х2, 0
(0= / |р ■ и Х_иХ21. (12)
—Х1, 0 —Х2, 0
Производство энергии возмущений ^(¿) рассчитывалось на основе полученного в [1, 9] интегрального уравнения энергетического баланса
где
XI, 0 Х2, 0
= — / с1:/ (1x2 Р и' и' -¡-^; </;
— XI 0 — Х2, 0
XI, 0 Х2, 0
/ р'
— XI, 0 — Х2, 0
(13)
ди XI ди Х2 5ж1 дж2
Л
XI, 0 Х2, 0
/ /
ди
XI
дж1
+
ди
XI
дж2
+
ди
X2
дж1
+
+
дм!
X2
дж2
ди' ди'
XI
X2
дж1 дж2
XI, 0
0
4 —
диди'
дж1 дж2
— XI, 0 —X2, 0
Здесь величины и X1, и X2, р' — пульсации компонент скорости и давления, ж1, о — п/в, ж2,0 — 20 Волновые чиела в для различных значений режимных параметров ^ьгь, М и т.„гЬ брались из табл. 1. Интегралы (г — 1, 2, 3, 4) в уравнении (13) описывают соответственно обмен энергией между возмущением и основным потоком, работу при пульсациоппом расширении или сжатии газа и вклады диссипативных процессов, определяемых сдвиговой и объемной вязкоетями. Так как коэффициенты переноса в системе (1) принимались постоянными, то при выбранной форме уравнения энергетического баланса вклад в него релаксации колебательной моды осуществляется косвенным образом через плотность и давление газа, связанные уравнением состояния с температурой.
Пульсацпонные характеристики течения определялись как разности
Ф '(¿,Ж1,Ж2) — Ф(г,Ж1,Ж2) - Ф8(г,Ж1,Ж2),
где компоненты вектор-функция Ф представляют собой мгновенные значения характеристик возмущенного течения, а Ф5 соответствующие характеристики несущего потока. Так как в данном случае в отличие от модельной задачи [9] невозмущенное течение не является точным стационарным решением системы (1), его мгновенные характеристики рассчитывались параллельно с расчетом возмущенного потока.
Интегралы вычислялись по квадратурным формулам трапеций с шагом к — 0.1 па сетке, использованной в расчетах по схеме (10). Для контроля вычислений значений производства энергии из энергетического уравнения (13) величины ^(¿) рассчитывались также непосредственно с использованием (11) и симметричной конечно-разностной аппроксимации временной производной:
Е(г + т)-Е(г-т)
27 '
Как показали расчеты, значения ^(¿), вычисленные двумя способами, различались не более чем на 1.5%.
2
2
2
-XI 0
— X2. 0
Примеры временных зависимостей кинетической энергии возмущений E(t, yvib) для некоторых режимов представлены на рис, 4, Видно, что эволюция E(t, a1; 7vib) соответствует кинематике развития структуры, отмеченной в комментарии к рис, 2, Из графиков следует: чем больше значение параметра 7vib, тем меньше кинетическая энергия структуры на всем временном интервале. Возрастание глубины возбуждения колебательной моды, как и увеличение объемной вязкости, приводит к дополнительной диссипации кинетической энергии возмущений. Временные зависимости модуля рейнольд-совых напряжений |axix2 |(t, Yvib, ai), определенных в (12), практически повторяют приведенные графики E(t, Yvib, a1). Расчеты также показали, что варьирование значения параметра Tvib от нуля до единицы слабо влияет на поведение временных зависимостей E(t, Yvib, ai) и |axix2|(t, Yvib, ai)-
Расчетные зависимости производства пульсациоппой энергии D(t), полученные из энергетического уравнения (13), приведены на рис, 5, Сопоставление графиков на рис, 4, 5 показывает, что на временном интервале, где происходит рост структуры и кинетической энергии возмущения, ее производство сначала положительно и возрастает, достигая некоторого максимума, а затем начинает убывать и в итоге становится отрицательным, что соответствует убыванию энергии структуры. Точки перехода кривых D(t) через нуль соответствуют максимумам па кривых E (t), а точки максимума и
D(t) E(t)
Для количественного сравнения вклада колебательной релаксации в диссипацию кинетической энергии возмущений находились средние по времени относительные отклонения
£в (ai,Yvib) = ©
© -1 /
E(t,a 1, Yvib) ~ E(t,ai, 0) E(t,a 1, 0)
■ 100% dt, Yvib = 0.1, 0.2, 0.4.
Рис. 4. Временные зависимости энергии возмущений Е(Ь, Yviь) для И,е = 100, М = 0.5 и Тгаб = 1 (сплошная линия — а\ = 0, штриховая — а = 2; 1 — ЪгЬ = 0, 2 — ЪгЬ = 0.2, 3 — 1ьИЬ = 0.4)
Рис. 5. Временные зависимости производства пульсационной энергии D(t, Yvib) Для Re = 100, M = 0.5 и Tvib = 1 (сплошная линия — а1 = 0, штриховая — а1 = 2; 1 — Yvib = 0, 2 — Yvib = 0.2, 3 — Yvib = 0.4)
Таблица 3. Относительные отклонения eß (7mb), % Для Re = 100, Tvib = 1
м а\ = 0 (xi = 2
ЪгЬ = 0.1 ЪгЬ = 0.2 ЪгЬ = 0.4 ЪгЬ = 0.1 ЪгЬ = 0.2 ЪгЬ = 0.4
0.2 0.031 0.183 0.853 0.033 0.187 0.860
0.5 1.980 4.055 11.790 2.024 4.153 12.109
Осреднение проводилось по условному времени "жизни" структуры в = 6. Результаты расчетов относительных отклонений (7vib) для некоторых наборов параметров представлены в табл. 3, Из таблицы следует, что рост параметра Yvib при фиксированных значениях числа Маха M, параметра объемной вязкоети а и времени колебательной релаксации Tvib приводит к большей диссипации кинетической энергии возмущения, В частности, в отсутствие объемной вязкости при а = 0 среднее по времени относительное уменьшение энергии E(t, Yvib) для Yvib = 0.4 и M = 0.5 достигает = 11.79%, Полученное значение хорошо согласуется с данными работы [9] для диссипации энергии возмущений E(t, Yvib) за счет только релаксационного процесса колебательной моды при отсутствии в течении других диссипативных процессов, хотя в [9] рассматривалась модельная задача, в которой воспроизводилось лишь затухание структуры, что в данном случае соответствует убывающим ветвям кривых на рис, 4, Это подтверждает сделанный в [9] вывод о том, что достижимый вклад в диссипацию возмущений релаксации колебательной моды соизмерим с вкладом объемной вязкости, по крайней мере, в относительных величинах. Отметим также (см, табл. 3) однонаправленность воздействия всех трех факторов — сжимаемости (числа Маха M), объемной вязкости (параметра о^) и возбуждения колебательных степеней свободы (параметра Yvib), возрастание которых вызывает увеличение диссипации энергии возмущений.
Анализ эволюции отдельных слагаемых в правой части энергетического уравнения (10) показал, что их величины и изменения во времени полностью аналогичны соответствующим данным работы [1]. Так, поведение производства энергии возмущений D(t) почти целиком определяется вкладом интеграла J1 (t), характеризующим процесс обмена энергией между возмущениями и несущим потоком в уравнении (10), При этом, как и в [1], затухание кинетической энергии возмущения связано в первую очередь со сменой знака рейнольдеовых напряжений в подынтегральном выражении J1(t), а не с отрицательными вкладами J3(t), J4(t) на всем времени "жизни" структуры. Следовательно, и в данном случае расчетное число Рейнольдса Re = 100 значительно превышает его критическое значение, которое может быть найдено из решения вариационной задачи энергетической теории устойчивости, вытекающей из уравнения (10) при D(t) = 0, Подобный вывод, сделанный в работе [1], нашел количественное подтверждение в расчетах вариационной задачи в работе [10], где были получены значения Recr < 10,
Из графиков на рис, 4, 5 и табл. 3 следует, что абсолютная величина расслоения данных при изменении параметра Yvib в предел ах 0 ^ 0.4 существенно меньше, чем при переходе от a1 = 0 к a1 = 2, Вместе с тем, как известно из кинетической теории газов [2], объемная вязкость (коэффициент a1 = 0) есть проявление слабо неравновесного процесса релаксации вращательных степеней свободы газовых молекул, В то же время значение параметра Yvib = 0.4 соответствует глубокому возбуждению колебательных степеней свободы [9], релаксация которых, как представляется, должна бы привести к более значительному диссипативному эффекту, чем объемная вязкость. Можно предположить, что выявленное несоответствие объясняется принятым в расчетах допуще-
нием, что диееипативные коэффициенты в системе (1), прежде всего вязкостные, не зависят от статической температуры потока, В результате сброс больших квантов колебательной энергии на поступательные степени свободы молекул, который вызывает повышение температуры потока, лишь косвенно, через связанные с температурой термодинамические параметры, учитывается в уравнении (10), Сделанное предположение будет проверено в дальнейшем в вычислениях с учетом температурной зависимости диссипативпых коэффициентов.
Заключение
В работе выполнено численное моделирование нелинейного развития дозвуковых вихревых возмущений в эволюционирующем во времени сдвиговом слое колебательно неравновесного газа. Диапазоны параметров течения, в частности, отношения коэффициентов объемной и сдвиговой вяз костей а и степени неравновесности колебательной моды Yvib соответствовали реальным значениям для двухатомных газов. Полученные результаты позволяют сделать следующие выводы,
1, Влияние колебательной релаксации на кинематику вихревой структуры с преимущественно невязким механизмом эволюции относительно мало, но ее дополнительный диссипативпый эффект проявляется вполне определенно,
2, При возрастании степени неравновесности колебательной энергии молекул газа в реальных для двухатомных газов пределах в отсутствие объемной вязкости относительное увеличение диссипации кинетической энергии структуры, осредненное по времени ее "жизни", достигает примерно 12%, что хорошо согласуется с аналогичными результатами исследований на модельной задаче [2].
3, Слабый по сравнению с влиянием объемной вязкости дополнительный диссипативпый эффект релаксации даже при глубоком возбуждении колебательной моды можно объяснить тем, что в использованной модели двухтемпературной гидродинамики не учитывалась зависимость диссипативпых коэффициентов, прежде всего вязкостных, от статической температуры потока.
Список литературы
[1] Григорьев Ю.Н., Ершов И.В., Зырянов К.В. Численное моделирование волн Кель-вина—Гельмгольца в слабо неравновесном молекулярном газе // Вычисл. технологии 2008. Т. 43, № 5. С. 25-40.
[2] Жданов В.М., А л невский М.Е. Процессы переноса и релаксации в молекулярных газах. М.: Наука, 1989.
[3] Осипов А.И., Уваров A.B. Неравновесный газ: Проблемы устойчивости // Успехи физ. наук. 1996. Т. 166, № 6. С. 639-650.
[4] Blumen W. Shear laver instability of an inviscid compressible fluid //J. Fluid Mech. 1970. Vol. 40, pt 4. P. 215-239.
[5] Michalke A. On the inviscid instability of the hyperbolic-tangent velocity profile // Ibid. 1964. Vol. 19. P. 543-556.
[6] Ковеня В.M., Яненко H.H. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, Сиб. отд-ние, 1981.
[7] Григорьев Ю.Н., Ершов И.В., Зырянов К.И., Синяя A.B. Численное моделирование эффекта объемной вязкости на последовательности вложенных сеток // Вычисл. технологии. 2006. Т. 11, № 3. С. 36-49.
[8] Квасов Б.И. Интерполяция кубическими и бикубическими сплайнами: Уч. пособие. Новосибирский гос. ун-т, 2004.
[9] Григорьев Ю.Н., Ершов И.В., Ершова Е. Е. Влияние колебательной релаксации на пульсационную активность в течениях возбужденного двухатомного газа // Журн. прикл. механики и техн. физ. 2004. Т. 45, № 3. С. 15-23.
[10] Григорьев Ю.Н., Ершов И.В. Энергетическая оценка критических чисел Рейнольдса в сжимаемом течении Куэтта. Влияние объемной вязкости // Там же. 2010. Т. 51, № 5. С. 59-67.
Поступила в редакцию 18 мая 2010 г., с доработки — 31 мая 2010 г.
