Stability of the Couette flow of a vibrationally nonequilibrium of molecular gas. Energy approach.pdf В работе [1] выведено уравнение энергетического баланса для плоскопараллельных течений колебательно-возбужденного молекулярного газа. На основе этого уравнения получено асимптотическое решение вариационной задачи определения критического числа Рейнольдса ламинарно-турбулентного перехода Rec в течении Куэтта колебательно-возбужденного молекулярного газа. В пределе малых волновых чисел возмущений найдена зависимость главного члена асимптотики от числа Маха M, коэффициента объемной вязкостей Пь, степени неравновесности колебательной энергии Yvib и времени колебательной релаксации tVt в виде Rec «1 + 3 где а1 = пЬ/П, П - коэффициент сдвиговой вязкости, Pr - число Прандтля, y - показатель адиабаты. 20 Y*bPr>/«1 + 4/3 (Y -1) M2 Pr' 1 - 1 - 10 ytVt + 33PrJ a + 4/3 2 п2 В данной работе вариационная задача, сформулированная в [1], решается численно во всем диапазоне волновых чисел продольных и поперечных возмущений. Течение Куэтта колебательно-возбужденного молекулярного газа описывается системой уравнений двухтемпературной аэродинамики, в которых учитывается зависимость коэффициентов переноса от температуры потока. В качестве температурной зависимости коэффициентов переноса выбран степенной закон Tn с показателем n < 1. Выбранная зависимость соответствует условиям относительно холодного несущего потока (мягким потенциалам межмолекулярного взаимодействия [2-4]). Предполагается, что удельные теплоемкости не зависят от статической и колебательной температур потока и постоянны. В соответствии с физическими представлениями [3, 4] модель двухтемпературной аэродинамики является общепринятой физико-математической моделью течений колебательно-возбужденного молекулярного газа, когда диссоциацией молекул, возбуждением верхних колебательных уровней и поправками на ангармонизм колебаний можно пренебречь. Основные уравнения и энергетический функционал Задача устойчивости течения Куэтта колебательно-возбужденного молекулярного газа рассматривается в расчетной области Q, представляющей собой прямоугольный параллелепипед, грани которого параллельны координатным плоскостям декартовой системы (xb x2, x3), а центр совпадает с началом координат. Непроницаемые бесконечные пластины, вдоль которых направлено основное течение, перпендикулярны оси x2. В качестве характерных величин для обезразмеривания использованы полуширина канала L по оси x2, модуль скорости потока U0 на непроницаемых стенках канала, постоянные плотность р0 и температура T0 основного течения и образованные из них время р0 = L/U0 и давление p0 = p0U02. Коэффициенты переноса обезразмеривались на их значения при температуре T0: сдвиговая и объемная вязкости на По и nb,0, а коэффициенты теплопроводности, обусловленные упругими энергообменами между поступательными степенями свободы молекул и неупругими обменами энергией вращательных и колебательных степеней свободы молекул с поступательными модами молекул, соответственно на V,0, Xrot 0, X,ib,0. В безразмерных переменных система уравнений двухтемпературной аэродинамики с учетом зависимости коэффициентов переноса от температуры потока записывается в виде = 0, Ф + д Р u д t д Xi ( \ д u, 1Л д д uд u. -L + u ,-i v д t J д x, v J д u д x, д p 1 д + дx Re дx, П(Т) П(Т) а, +- Р 3 )д х, д x Re J J дТ дТ Л 1Ч д щ y PllT + u IT 1 + (Y-1)РТ= дТ д x, П(Т) дt дx J дxt RePr дТ дТ дТ Ya2 vib vib vib П(Т )д x. д x, RePr Yvib Р + u: (1 - Yvib) V дt + Yvib Р ^vib -Т) tvt(1 - Yvib) Yvib Р (ТУ!Ъ - Т) TVT(1 - Yvib) (1) YM2p = РТ, п(т)= Тп, i = 1,2,3, j = 1,2,3, где p, u, p, Т, тл, tvt - плотность, компоненты вектора скорости, давление, статическая и колебательные температуры газа и время колебательной релаксации соответственно, а по повторяющимся индексам подразумевается суммирование. В уравнении энергии системы (1) опущена группа нелинейных слагаемых, составляющих диссипативную функцию. Такое приближение является распространенным в задачах устойчивости сжимаемых течений [5]. Параметры, входящие в уравнения системы (1), определяются следующим образом: ц = nb0/ По - отношение объемной и сдвиговой вязкостей; а2 = 1 vib,0 /(Хtr, 0 + 1 rot, 0 ) ; параметр Yvib = cvib /(ctr + Crot + Cvib ) определяет долю внутренней энергии газа, приходящуюся на колебательные моды молекул, и в каком-то смысле характеризует степень неравновесности последних [6-8]; безразмерные критерии Re = U0LP0/П0 , M = U0NYRT0 и Pr = П0 (ctr +crot)/(xtr,0 + Xrot,0) есть соответственно числа Рейнольдса, Маха и Прандтля несущего потока, где Y = (ctr + crot + R)/(ctr + crot) - показатель адиабаты, R - газовая постоянная и ctr, crot, cvib - соответственно удельные теплоемкости при постоянном объеме, определяющие энергоемкость поступательных, вращательных и колебательной мод молекул газа. Система (1) описывает распространенную в аэродинамике ситуацию, когда характерные времена микроскопических процессов энергообмена между различными степенями свободы молекул газа оцениваются системой неравенств [3, 4]: ttt ~ trt • • • > "n ' "v,^ uv,j ' • • •' uv, N r (x2) = (м0, M, w, w,. -, w VN N а матрицы G, F размером 5(N+1)x5(N+1) вычисляются с использованием специальной процедуры Matlab по формулам G = A ® DN + B ® DlN + C ® /N, F = K ® /N, где знак ® обозначает прямое (тензорное) произведение матриц [11]; IN - единичная матрица размером (N+1)x(N+1); A, B, C - матрицы размером 5x5: A = Г1 0 0 0 0 0 4 а1 + 0 0 0 13 0 0 1 0 0 0 0 0 1 0 0 0 0 0 20 y vib 33(1-Yvib). B = 0 -(а2 + S2) C = -ал (Y-1)M2Pr ( -|52 +(aj +4 |а2 0 -aSl а, +1 0 /а| а1 + з - 2(y-1)M2 Pr 0 V - а5| а, + з 0 -| а2 + 4]б2 0 iS| а, +0 0 0 /ал ~2~ 0 -(а2 + S2) 0 П о 0 /5|а, +3I 0 0 0 0 10 0 0 20yль (а2 +52) 33(1- Yvib) ( 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 Yvib Pr - Yvib Pr Y(1 - Y vib) т vt Y vib Pr Y(1- Yvib) tvt Yvib Pr 0 0 0 - Y(1-Yvib) T VT Y(1- Y vib) T VT Однородные граничные условия (7) для уравнения (10) учитываются неявно через оператор DDN и на дискретном уровне реализуются заменой матриц DN (к = 1, 2) на окаймленные матрицы размером (N-1)x(N-1) [9, 10]. Последние получаются при выполнении условий D0,j = DN,j = 0, DI0 = dIn = 0, i = 0,1,.,N, j = 0,1,.,N, I = 1,2. Для нахождения всех собственных значений и функций обобщенной спектральной задачи (10) использовалась процедура Matlab, реализующая QZ-алго-ритм [12], который позволяет одновременным ортогональным преобразованием привести пару матриц G, F к обобщенной верхней треугольной форме. K = В результате применения данной процедуры для фиксированных значений числа Маха M, объемной вязкости а1, степени неравновесности колебательной энергии у т;ь, времени колебательной релаксации tVt и каждой пары волновых чисел (а, 5) получался набор Л/+1 собственных значений, среди которых находилось минимальное по модулю число Рейнольдса Re(a,5) =2 | Xmin(a,5)|. Значение критического числа Рейнольдса Rec для данных M, а1, tVt и у vib принималось равным минимальному значению Re во всем диапазоне волновых чисел Re(a,5): Rec = min Re(a, 5). Затем вычислялись соответствующие Rec собственные функ- (a 5) ции u, v, w, 9, 9v. Вычисления спектров собственных значений Х(а, 5, M, а1, tVt, yvib) выполнялись для случая, когда в качестве молекул несущего газа рассматриваются «мак-свелловские» молекулы, тогда n =1 [2, 3]. Все расчеты велись в диапазоне волновых чисел а = -5^5, 5 = -5^5 при следующих значениях параметров: yvlb = 0,1^0,4; tVt = 1^4; ai = 0^2; M = 0,1^1; Pr = 3/4; у = 7/5. Шаги изменения волновых чисел были выбраны равными ha = h5 = 0,1. Число узлов коллокации в интервале х2е[-1, 1] принималось равным Ж+1 = 50. Для проверки точности расчетов проводилось варьирование числа узлов коллокации в диапазоне Л/+1 = 32^100. Результаты расчетов и их обсуждение Расчеты показали, что при всех рассматриваемых значениях степени неравновесности колебательной энергии yvib, времени колебательной релаксации tVt, объемной вязкости а1 и числа Маха M минимальные по модулю собственные значения Re(a,5) =2 | Xmin(a,5)| достигаются на оси а^0 (при 5=0) в плоскости волновых чисел (а, 5). Изолинии Re(a,5) приведены на рис. 1. 5 5 Рис. 1. Линии уровня поверхностей Re(a,5) для M=0,5, a1=0, Tvr=1 (а - Yvib=0,2; б - Yvib=0,4) Как и в случае умеренного возбуждения, рассмотренном в [13], наиболее «опасными» являются возмущения продольной моды. С учетом периодичности полученного решения по продольной координате х1 эти возмущения представляют собой пары двумерных вихрей, вращающихся в противоположных направлениях, с осями, перпендикулярными несущему потоку. Распределение завихренности в этих вихрях вычисляется по формуле ra(xj,x2) = -| avi + ^ЬЛcos(aXj)-| avr -sin(axj). I dx2 ) \ dx2 ) Здесь ur (x2), ui(x2), vr(x2), vi(x2) - вещественные и мнимые части собственных функций u, v. На рис. 2 представлены изолинии завихренностей ra(xb x2), при различных критических числах Рейнольдса Rec(a, M, ab tVt, yvib) и значениях амплитуд возмущений скорости, составляющих 10 % значения модуля скорости несущего потока на непроницаемых границах. Рис. 2. Изолинии завихренности ro(x1, x2) для M=0,5, a1=0, tVt=1 (а - yvib=0,2, Rec=8,56; б - yvib=0,4, Rec=9,66) Re Re Рис. 3. Зависимости числа Рейнольдса Re(a) для продольных мод возмущений при M = 0,5 (а - a1=0; б - a1=2; 1,1' - yvib=0,2; 2, 2' - yvib=0,3; 3, 3' - yvib=0,4; штрихпунктирная линия -зависимость критического числа Рейнольса Rec от волнового числа a) Зависимость числа Рейнольдса для продольных мод возмущений от волнового числа a представлена на рис. 3, где штрихпунктирные линии соединяют значения абсолютных минимумов на параметризованных по yvib и tVt кривых Re(a), что позволяет проследить эволюцию Rec. На рис. 4 приведена зависимость Rec от степени неравновесности yvib. На рис. 1, 3, 4 видно, что с увеличением значений параметров M, a1, tVt, yvib критические числа Рейнольдса Rec и соответствующие им значения волнового числа а возрастают. Таким образом, результаты расчетов показывают, что представленные выше асимптотические выражения (9) качественно описывают зависимость Rec от параметров течения. Критические значения числа Рейнольдса Rec(a1,тvT,yVlb,M) приведены в табл. 1, а соответствующие им значения волнового числа а в табл. 2. В рассмотренном диапазоне значений параметров задачи число Рейнольдса Rec в пределе увеличивается примерно в два с половиной раза. При этом влияние каждого параметра на Rec при фиксированных значениях остальных параметров существенно различается. Следует отметить, что наиболее значительное влияние на Rec оказывает степень неравновесности колебательной моды yvib. Оценки, выполненные в [6], показывают, что возбуждение колебательной моды yvib в рассмотренном диапазоне достаточно легко достигается с помощью лазера. Поэтому лазерная накачка колебательных мод может стать эффективным способом управления течениями молекулярных газов. Рис. 4. Зависимости критических чисел Рейнольса Rec от степени неравновесности колебательной моды yvib (а - а1=0; б - а1=2; 1, 1' - M=0,2; 2' - M=0,4; 3, 3' - M=0,6; 4, 4' -M=0,8) Таблица 1 Критические значения числа Рейнольдса Rec( M, а!, Тут, Yvib) M TVT=1 Tvt=4 Yvib=0,1 I Yvib=0,2 1 Yvib=0,4 Yvib=0,1 I Yvib=0,2 1 Yvib=0,4 а1=0 0,2 Rec=7,469 8,213 9,264 9,264 10,75 12,85 0,5 Rec=7,785 8,560 9,656 9,656 11,21 13,40 0,8 Rec=8,366 9,199 10,38 10,38 12,04 14,40 а1=2 0,2 Rec=9,614 10,57 11,92 11,92 13,84 16,54 0,5 Rec=10,02 11,02 12,43 12,43 14,42 17,24 0,8 Rec=10,77 11,84 13,36 13,36 15,50 18,53 Таблица 2 M TVT=1 Tvt=4 Yvib=0,1 1 Yvib=0,2 1 Yvib=0,4 Yvib=0,1 1 Yvib=0,2 1 Yvib=0,4 а1=0 0,2 а=0,482 0,526 0,613 0,613 0,788 1,139 0,5 а=0,488 0,533 0,621 0,621 0,799 1,154 0,8 а=0,554 0,604 0,705 0,705 0,906 1,309 а1=2 0,2 а=1,111 1,211 1,413 1,413 1,817 2,625 0,5 а=1,125 1,228 1,432 1,432 1,842 2,660 0,8 а=1,277 1,393 1,625 1,625 2,089 3,017

Скачать электронную версию публикации