Численное моделирование эффекта объемной вязкости на последовательности вложенных сеток Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 11, № 3, 2006

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЭФФЕКТА ОБЪЕМНОЙ ВЯЗКОСТИ НА ПОСЛЕДОВАТЕЛЬНОСТИ ВЛОЖЕННЫХ СЕТОК*

Ю. Н. Григорьев Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: [email protected]

И. В. Ершов, К. И. Зырянов

Новосибирская государственная академия водного транспорта, Россия

e-mail: [email protected]

А. В. Синяя

Новосибирский государственный университет, Россия

The influence of bulk viscosity on evolution of a finite-amplitude vortex disturbance riding on a shear flow of molecular gas is studied in the framework of full Navier — Stokes equations. The result of calculations on a sequence of embedded meshes have shown that the effect of scheme viscosity in energy dissipation is negligible compared to the bulk viscosity.

Введение

Диссипативный эффект, возникающий в молекулярных газах в процессе термической релаксации, известен с середины 30-х годов прошлого столетия [1, 2]. В последнее время он привлек внимание аэродинамиков в связи с перспективой использовать его для затягивания ламинарно-турбулентного перехода и подавления турбулентности. Как известно [3], при умеренном термическом возбуждении этот эффект описывается коэффициентом объемной вязкости в уравнениях Навье — Стокса.

В [4] проведены сравнительные эксперименты по ламинарно-турбулентному переходу в течении Пуазейля в круглой трубе на азоте N2 и окиси углерода CO. Эти газы почти идентичны по физическим свойствам, но объемная вязкость CO, рассчитанная по данным о затухании ультразвука, в несколько раз превышает аналогично вычисленную величину для N2. В экспериментах установлено, что при одинаковых условиях число Рейнольдса перехода в более "вязком" CO примерно на 10 % превышает соответствующее значение для N2. По ряду причин достоверность полученных данных представляется дискуссионной.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-01-00359).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

В работе [5], отчасти инспирированной публикацией [4], в рамках линейной теории устойчивости численно исследовано влияние возбуждения внутренних степеней свободы молекул на ламинарно-турбулентный переход в сжимаемом пограничном слое на пластине. Однако в расчетах для сверхзвукового воздушного потока получено, что изменение числа Рейнольдса перехода при учете объемной вязкости не превышает десятых долей процента. Аналогичный результат дали наши оценки [6], полученные для предельно допустимого в рамках модели Навье — Стокса диапазона чисел Маха и отношения коэффициентов объемной и динамической вязкостей.

Тем не менее эти результаты, полученные в линейном приближении, не могли служить прямым опровержением экспериментов [4], поскольку в них переход к турбулентности исследовался вплоть до нелинейной стадии, тогда так в линейном приближении течение Гагена — Пуазейля устойчиво.

Для рассмотрения нелинейного этапа перехода в [7] предложено исследовать модельное течение с начальными данными, составленными суперпозицией вихря Рэнкина с линейным сдвиговым потоком. Несмотря на простоту, такая модель вполне адекватно воспроизводит характерный момент современных сценариев ламинарно-турбулентного перехода и генерации развитой турбулентности, состоящий в эволюции крупных вихревых структур на фоне несущего сдвигового течения [8]. Проведенные в [7] расчеты такого течения на основе полных уравнений Навье — Стокса показали, что в реально достижимом диапазоне значений объемной вязкости диссипативный эффект, оцениваемый по скорости затухания начального вихревого возмущения, достаточно существен.

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

1. Постановка задачи и основные уравнения

В модельной ячейке, представляющей квадратную область в поле течения, начально-краевая задача ставилась следующим образом. Динамика вихревого возмущения описывалась системой полных уравнений Навье — Стокса сжимаемого теплопроводного газа. В качестве характерных величин для обезразмеривания выбраны начальный диаметр структуры 2До, время т0 = 2Д0/и0 и давление р0 = д0^о, а также модуль скорости и0, плотность д0 и температура Т0 на границах расчетной области, параллельных основному потоку. Завихренность обезразмеривалась на величину = 2 и0/1. В безразмерных переменных система уравнений для построения разностной схемы записывалась в операторном виде:

!+А=и (1>

Здесь ^ = (д,их,иу, Т) — вектор решения, компоненты которого есть плотность, скорости в направлении координатных осей х,у и температура соответственно. Операторы А и И,

в уравнении (1) имеют вид

(

А

а1

Т д 7дМ2 дх

Т д 7дМ2 ду

^ дх а1 + а2

0

йду 0

а1 + аз

1 д_ ^М2 дх

1 д ^М2 ду

дд 0 ^ - 1)ТдХ (7 - 1)Тду а1 + а4

(2)

д

д

а1

аз

их п + иу п ,

дх ду (4 +За) д2

а2

(4 +За) д2 1 д2

ЗдИе дх2 дИ,е ду2 1 д2 7 д2

ЗдИе ду2 дИ,е дх2

а4

1

д2

дИеРг дх2 дИеРг ду2

Я,

4 —

7(7 - 1)М2

дИ,е

Ях

дих

ду

И. — (0, Ях, Яу, Яа), (1 + За) д2иу (1 + За) д2их

(3)

ЗдИе дхду'

Яу

ЗдИе дхду'

+

диу

дх

+

дих

дх

+

диу

ду

+ |а- З

2 \ I дих диу

+

дх ду

Предполагалось, что диссипативные коэффициенты в системе (1) не зависят от температуры и постоянны. Градиенты давления в уравнениях импульсов были исключены с помощью уравнения состояния идеального газа, которое в безразмерных переменных имеет вид р — дТ/'уМ2.

Для системы (1) на границах расчетной области во все моменты времени выполнялись условия:

при х — ±х/2, у е [-х/2, х/2]

ихх/2, у) — их^, -х/2, у), иух/2, у) — -иу-х/2, у),

д(Ь, х/2, у) — д(Ь, -х/2, у), Т(г, х/2, у) — Т(Ь, -х/2, у); при у — ±х/2, х е [-х/2, х/2]

их(Ь, х, х/2) — -их(Ь, х, -х/2), иу (Ь, х, х/2) — иу (г, х, -х/2), д(Ь, х, х/2) — д(г, х, -х/2), Т(г, х, х/2) — Т(г, х, -х/2).

(4)

(5)

Несущий поток в модельной ячейке задавался как точное стационарное решение системы (1)-(3) (течение Куэтта) с граничными условиями (4), (5):

им — *ТМ — 1 + И - 1)М2РГ

х

2

1

4у2 х2

0

2

2

2

2

^s(y) = у е [ х/2, х/2].

(6)

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

пх(0, х, у)

Пу(0, х, у)

2у/х + ву/2(х2 + у2), х2 + у2 > 1/4, 2у/х + 2ву, х2 + у2 < 1/4,

-вх/2(х2 + у2), х2 + у2 > 1/4, -2вх, х2 + у2 < 1/4,

(7)

а начальное распределение термодинамических величин соответствовало условиям невозмущенного потока (6).

Параметры, входящие в уравнения (1)-(3) и условия (4)-(7), определяются следующим образом. Коэффициент а равен отношению объемной вязкости к сдвиговой а = ^ он характеризует степень неравновесности внутренних степеней свободы молекул газа. М = Цо/л/тЯТо есть число Маха несущего потока, И,е = 2иоЯо— число Рейнольдса, Рг = ^Ср/А — число Прандтля. Здесь Я есть газовая постоянная, 7 = ср/с„ — показатель адиабаты, ср и а€ — удельные теплоемкости при постоянных давлении и объеме соответственно, А - коэффициент теплопроводности. Параметр в = ^оЯо/(2Цо) определяет относительную интенсивность вихревого возмущения. Величина х = 1/(2Яо) характеризует коэффициент перемежаемости в возмущенном потоке.

Все расчеты проводились при следующих значениях параметров: М = 0.5; И,е = 100; Рг = 0.74; в = 0.2; х = 3; 7 = 1.4; а = 0 ... 2.

2. Разностная схема и метод решения

В расчетах использовались регулярные сетки с одинаковым шагом к по обеим пространственным координатам. Система (1) аппроксимировалась весовой конечно-разностной схемой с расщеплением по физическим процессам и пространственным переменным [2]. В операторной форме схема имеет вид

сга+1 _ си

--Н + л!2) + (1 - ОД] = К. (8)

т

Здесь СП = ($, «X н, «П н' Ти,) — сеточная вектор-функция решения на п-м временном

(2)

слое в узле (г, т — шаг по времени; 5 — весовой параметр. Оператор АН ) включает симметричные аппроксимации со вторым порядком первых и вторых пространственных производных на трехточечном шаблоне по каждой пространственной координате и записывается следующим образом:

А(2) = 0(2) + 0(2) + А(2) АН = их, Н + ^у,Н + Ах ,Н,

где Б^Ь, I — х, у — диагональные оператор-матрицы вида

I и7,НЛ£

(2)

Б

(2)

0 0

ип л(2) л

Щ.ь+ч - а1

ипл? - а2

а1

л

2)

11

ЗИед;

а оператор-матрицы (

А

2)

Тпя Т и и.

0 0

[(4 + За)5х,г + З5уЛ], Л2

0 дП л[р5х,е

0

и

л 2)

1

ч, ь1 ^

-л(2)

ИеРгдп"11

л 2) л11

ЗИед

[З5х,1 + (4 + За)иу,1],

дП л?\г

ь °х,£ л ( 2) 2 л1

1йП М

Ткиу,е_ л(2)

их,£ л(2) 1 ле

7М

21 Ч

1дП М

21 Ч

V

0

(7 - 1)ТЬл12)их,1 (7 - 1)ТПл12)5у,е

5у,1

7М

л(2)

21 Ч

0

Здесь 5Ы ются как

символ Кронекера; симметричные разностные операторы л^2) и лЦ) выража-

л(1) +л(1) л(2) — л1 (+) +л1 (-) , л(2)

2

л(1) л(1)

л1 (+) л1 (-)

2

Индексы " +" и "-" означают направленные разностные отношения в узле (г, ]) "вперед" и "назад" соответственно.

Оператор ИП в правой части равенства (8) состоит из симметричных по каждой пространственной координате аппроксимаций со вторым порядком смешанных производных из уравнений импульсов и слагаемых диссипативной функции из уравнения энергии и записывается в виде

т>га _ (п тэ п тэ п тэ п \

ИЬ — (0, Ях,Ь, Яу,Ь, Я4,Ь),

ях ь — (^За1лх2)лу2)ип яп ь — (1+За)лу2)лх2)

7(7 - 1)М2

ЗИедп иу'ь ±Ьу'н

ЗИедп "у х х'ь

тэп

И,ед^

(лу2)их,ь + лх2)ип,ь)2 +

у,Ь) 2

+2(лх2)их,ь)2 + 2(лу2)ип,ь)2 + [а - ^ | (л^ + л^н)2

Определенная таким образом разностная схема (8) аппроксимирует исходную дифференциальную систему (1) с порядком 0(т + к2) и абсолютно устойчива при 5 > 0, 5.

0

0

0

0

0

0

0

0

0

0

0

П

Следуя [2], после записи схемы (8) в каноническом виде и приближенной факторизации

оператора ^Е + совершается переход к эквивалентной разностной схеме того же

порядка аппроксимации, что и (8), которая реализуется дробными шагами в следующем виде:

Е + тЖ^) 1 = -А^Г + Ип, Е + тЖ^ )С+1 = 1,

Е + т^)^4 = £п+ 2, (9)

Е + т^) С+1 = С+3,

£«,+ 1 = £П + т £«+1

сп+т (с«+т с«+т сп+т сп+т) , 0 ч -

где £п+ 4 = 4 , 4 , 4 , 4 I, т =1, 2, 3, есть вектор приращений сеточных

вектор-функций на промежуточных шагах (п + т/4) относительно их значений на слое п, а Е — единичная матрица. Алгоритм реализации схемы в дробных шагах описан в [10].

3. Тестирование схемы

Эффект дополнительного подавления возмущений, связанный с объемной вязкостью, предварительно оценивался на уровне нескольких процентов. Поэтому было необходимо, чтобы вычислительные погрешности не превышали величин третьего порядка малости. С этой целью использованная в расчетах разностная схема (9) тщательно тестировалась. Для тестирования взяты точные решения системы Навье — Стокса (1), которые являются составными элементами рассматриваемого модельного течения.

Параметры течения в тестовых расчетах были идентичны использованным в основной серии. Шаг по времени во всех расчетах т = 0.01. Рассматривалась последовательность вложенных сеток с шагом h = 0.1; 0.05; 0.025; 0.0125. Наиболее мелкая сетка содержала 241 х 241 узлов. Это позволило определить порядок практической сходимости реализованной в работе схемы.

Все расчеты выполнялись с удвоенной точностью. Относительные погрешности оценивались в нормах пространств сеточных функций L2 , h и Ch по формулам

\\Ph< \\L2,h \\Ph< \\Ch

i V/2

где \\<h\U2 h = E I<ij |2 , \\<h\\ch = max ||; оператор проектирования на сетку Ph

\(i ,j) J (i ' j)

на точных решениях определялся как Ph<(x, y, t) = <(ih, jy, пт).

3.1. Стационарный несущий поток. Течение Куэтта

Возмущения, временная эволюция которых прослеживалась в модельном течении, определялись следующим образом:

<h = $h - Ph$0, (11)

Рис. 1. Временные зависимости погрешностей ||ДиУ)н||сй: ◦ — Н = 0.1; * — Н = 0.05; +--Н = 0.025;

• - Н = 0.0125.

Таблица 1. Установившиеся значения погрешностей аппроксимации точных решений (6) при £ = 8 для различных значений шага сетки Н

Погрешность Н = 0.1 Н = 0.05 Н = 0.025 Н = 0.0125

Дпх,н • 106 8.5190 1.5210 0.3405 0.0951

— 5.6009 4.4670 3.5804

Д-у,н • 106 2.1660 0.4200 0.2580 0.1080

— 5.1571 1.6279 2.3889

Д^н • 106 6.5970 1.3780 0.3110 0.0866

— 4.7874 4.4309 3.5912

ДТн • 106 0.4513 0.0946 0.0213 0.0058

— 4.7706 4.4413 3.6724

где Фн — мгновенные значения расчетных сеточных характеристик; Фо — характеристики стационарного несущего потока (6). В этой связи было необходимо, чтобы разностная схема с высокой точностью сохраняла стационарные профили (6), являющиеся точным решением задачи Куэтта для системы Навье — Стокса (1) с граничными условиями (4), (5).

Расчеты проводились на временах Т = 800т. Зависимости от времени вычислительных погрешностей для пх>н, иу,н, Он и Тн проявляли универсальное поведение. Погрешности на наиболее грубой сетке 31 х 31 значительно превышают погрешности на более мелких сетках. На всех сетках отклонения от точного решения выходят на стационарные значения за времена порядка 50т. Характерные зависимости от времени погрешностей иу,н в равномерной норме Сн приведены на рис. 1. Следует отметить, что именно || ДмУ)н||сь подвержены наибольшим флуктуациям, так как в данном случае фиксируется отклонение от нуля.

Установившиеся значения относительных погрешностей для всех величин в норме Сн даны в табл 1. Там же даны отношения погрешностей на двух последовательных сетках.

Можно констатировать, что, за исключением Диу,н, все параметры сходятся к точному решению приблизительно с теоретическим порядком аппроксимации данной схемы О (к2). Одновременно результаты показывают, что на интервале времени £ < 800т относительные отклонения от точного решения (6) для всех сеток не превышают величины порядка 10-6.

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

3.2. Диффузия завихренности. Вихрь Рэнкина

В качестве теста на нестационарном решении рассчитывалось течение, названное нами изохорным вихрем Рэнкина. При этом диффузия завихренности, первоначально равномерно распределенной в круге радиуса r < 1/2, происходит при постоянной плотности q = const. На границе расчетной области ставились условия (4), (5), а начальные условия соответствовали (7) при отсутствии сдвигового течения. В силу вращательной симметрии задачи расходимость поля скорости равна нулю и уравнение неразрывности в системе (1) удовлетворяется тождественно. Из уравнений импульсов исчезают слагаемые с объемной вязкостью, и они сводятся к уравнениям вязкой несжимаемой жидкости. Переходя от них к уравнению для скалярной завихренности, получаем для него известное точное решение [3]. В наших безразмерных переменных оно имеет вид

1/2

exRe ( r2 ReW ( Í 2Re\ ZrRe — exP (--—I I exp (--— ) Jo '

o

( t) exRe ( r2Re\ f ( £2Re\ T (rRe \

= —ex4J explv(12)

где /0(С) — модифицированная функция Бесселя нулевого порядка; г = л/х2 + у2 — расстояние от центра вихря до точки расчетной области с координатами (х, у). В силу принятого предположения о независимости коэффициента динамической вязкости от температуры уравнение энергии может быть решено отдельно и в данной работе не рассматривалось.

В приведенной ниже модельной задаче можно выделить два характерных времени. Одно из них определяется как время существования крупного вихревого возмущения в сдвиговом потоке. Обычно оно составляет 4-6 времен собственных оборотов вихря относительно центра завихренности [8]. Этим определяется верхний предел времени моделирования эволюции структуры. Для вихря Рэнкина с постоянной завихренностью П0 угловая скорость вращения е = П0/2. При этом период обращения составляет £0 = 2п/е. Безразмерное время оборота в нашем случае £0/т0 = п/в = 15, 7. С другой стороны, время счета в постановке задачи с периодическими условиями естественно ограничено моментом времени ¿е, когда значение вихревого возмущения на границе достигает некоторой малой величины.

В наших расчетах безразмерная завихренность вихревого возмущения в начальный момент при г < 1/2 была ^о(0, г) = 2хв = 1-2. Из физических соображений относительная величина возмущения на границе, определяющая окончание счета, принималась на уровне

и(х/2, ¿е)М = е(х/2) = 10-2.

С использованием точного решения (12) можно получить предварительную оценку времени ¿е следующим образом. В [3] дано асимптотическое разложение решения (12) в виде

w(r, t)

-= exp

Wo

(r2 + a2)Re

4t

ж k

E 2t «*) ■ (13)

k=1

справедливое для внешней области при г > а. Здесь а =1/2 — начальный радиус вихря Рэнкина; минимальный радиус граничных точек г = х/2 = 3/2; /к (£) — модифицированная функция Бесселя к-го порядка. Для оценки входящей в (13) суммы заметим, что /к(С) < /о(С), к =1, 2, ... [11]. Предварительно приняв £ ~ 10, имеем

V П k Tk (-Re) < To (-Re) = 9.1281 = 4.564.

^\r) kV2t J - 0V2t /1 - (a/r) 2

k=i w '

ж

После чего получаем

+ ~

Ье -

103

10.24.

(14)

16(21п 10+ 1п 4.564)

Видно, что ¿е не превосходит времени одного оборота вихря. Время окончания счета ¿е по уровню е(х/2) = 10-2 было уточнено в процессе тестов. Значения относительной завихренности в, коэффициента перемежаемости х, числа Рейнольдса И,е, определяющие величину ¿е (14), выбраны в физически обоснованных пределах. Поэтому полученная оценка показывает, что в рамках данной модели можно проследить лишь сравнительно небольшой отрезок естественной эволюции крупного вихревого возмущения.

Интеграл в формуле точного решения (12) вычислялся с помощью квадратурной формулы трапеций с шагом к^ = 0.001. По рассчитываемому скоростному полю сеточная функция завихренности определялась с помощью симметричной конечно-разностной ап-

проксимации

га = Х

ши = 2

их,г^-1 их,г, иу,г-1, ] иу,г+1, ]

2к

2к

(15)

Относительные погрешности оценивались в нормах (10). Необходимо отметить, что в данном случае кроме погрешностей, присущих собственно схеме (9), имеются другие источники ошибок. Они связаны с численным дифференцированием (15) и приближенным интегрированием в (12). В этой связи здесь нельзя было ожидать сходимости с порядком О (к2). Поэтому в данной работе определялись прежде всего качественные характеристики эволюции осесимметричного вихревого возмущения в квадратной области. В частности, установлено, что допустимая величина относительного возмущения на границе е(х/2) = 10-2 на самой грубой сетке с к = 0.1 достигалась при ¿е = 600т = 6. Сравнение с оценкой (14) показывает, что здесь возникает дополнительный диссипативный эффект вычислительной природы. Это определило временной интервал дальнейших расчетов.

Картина изолиний сеточной завихренности на плоскости (х, у) показала, что на временном интервале до ¿е = 6 диффузия завихренности происходит с сохранением высокой степени вращательной симметрии. Пример сопоставления расчетных данных (15) с точным решением (12) приведен на рис. 2.

Рис. 2. Эволюция завихренности в сечении плоскостью х = 0 для сетки с шагом Н = 0.025: сплошные кривые — расчет по формуле (12); значки — численные значения, полученные по формуле (15), * — £ = 2; +--£ = 4.

2

3

1 2 3 4 5 6t

Рис. 3. Эволюция во времени погрешностей аппроксимации завихренности Аи^ в норме С^: 1 — Н = 0.1; 2 — Н = 0.05; 3 — Н = 0.025; 4 — Н = 0.0125.

На рис. 3 показана эволюция во времени относительного отклонения Aw>h в равномерной норме. Видно, что на интервале времени [0, te] погрешности возрастают примерно на порядок. Тем не менее даже на самой грубой сетке погрешность остается в приемлемых пределах ||Awh|| < 5 • 10-3. Погрешности в норме L2,h проявляют аналогичное поведение. На усредненных по времени погрешностях отмечалась сходимость с уменьшением шага h. При этом отношение погрешностей на двух последовательных сетках было в пределах 1.4... 1.5 (сравните с табл. 1).

В качестве дополнительного теста рассчитывалась эволюция неизохорного вихря Рэн-кина при q = const. В этом случае в силу вращательной симметрии течения зависимость от времени завихренности в центре вихря одинакова для сжимаемой и несжимаемой жидкостей. Соответствующее аналитическое выражение [3] при нашем обезразмеривании имеет вид

w(0, t) = 2вх (1 - exp (-Re/16t)). (16)

Сопоставление численных результатов с аналитической зависимостью (16) свободно от большинства дополнительных погрешностей, возникающих при тестировании на точном решении (12). Результаты расчетов показывают, что на большей части временного интервала [0, te] на последовательности сеток имеет место сходимость к точному решению (16) с порядком O(h). Пример соответствующих данных приведен в табл. 2.

Таким образом, можно заключить, что снижение порядка сходимости в тестах с вихрем Рэнкина до O(h) связано с численным дифференцированием (15) сеточного поля скорости. Этот нежелательный эффект не имеет значения для расчетов модельного течения (4)-(7), в которых использовались только скоростные переменные.

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

Таблица 2. Относительное отклонение от точного решения (16) при t = 6

h = 0.1 h = 0.05 h = 0.025 h = 0.0125

||Д^|| • 103 2.784 1.722 0.999 0.569

— 1.617 1.724 1.756

4. Расчет модельного течения на последовательности вложенных сеток

Для оценки влияния объемной вязкости на пульсационные характеристики модельного течения, как ив [7], рассматривалась эволюция во времени кинетической энергии возмущения

х/2 х/2

2 у у (Ч ■ -у.

-Х/2 -х/2

Е (¿) = х / ¿ж / (<2 + <) (17)

и абсолютной величины рейнольдсовых напряжений

Х/2 х/2

(¿) = J ¿ж J ¿у (18)

-Х/2 -Х/2

Кроме того, рассчитывалась величина производства энергии возмущений на основе выведенного в [7] интегрального уравнения баланса

¿Е 1

Я(Г) = — = 71 + /2 - — (7а + ), (19)

аг Ке

где слагаемое

7 = -I -, ^

-Х/2 -Х/2

описывает обмен энергией между возмущением и основным потоком. Интеграл

Х/2 Х/2

I I , / дм' дмУ

= ¿ж / аур' + У

дж ду

Х/2 Х/2

интерпретируется как работа при пульсационном расширении (сжатии) газа. Интегралы

Х/2 Х/2

7а = J ¿ж J ¿у

-Х/2 - Х/2

д<л2 + (дмХЛ2 + ( д^УЛ2 + ( д^УЛ2 + !(дмХ + д<

дж у \ ду у \ дж у \ ду у 3 \ дж ду

Х/2 Х/2

Г с /дм' дмУ 74 = ¿ж ¿у +

дж ду

-Х/2 -Х/2

определяют вклады диссипативных процессов, связанных с динамической и объемной вяз-костями соответственно. В уравнениях (17)-(19) и', , р' — пульсации скорости и давления, определяемые по формуле (11); интегралы вычислялись по квадратурным формулам трапеций с шагом Н = 0.1; 0.05; 0.025. Расчеты выполнялись на последовательности вложенных сеток с шагами Н = 0.1; 0.05; 0.025.

2

2

На рис. 4 и 5 приведены зависимости от времени величин Е(Ь) и О(Ь) для сетки с шагом к = 0.025. Точки на рис. 5 получены на основе конечно-разностной аппроксимации

П(0 = ^ = Е(' + Т) - Е(' - Т) (20)

аЬ 2 т

с использованием рассчитанных значений (17). Как следует из приведенных зависимостей, с ростом коэффициента объемной вязкости а кинетическая энергия возмущений затухает все более интенсивно. Соответственно скорость диссипации кинетической энергии также возрастает. Для всех использованных сеток эти графики практически неразличимы и совпадают с аналогичными зависимостями, рассчитанными в [7] на сетке с шагом к = 0.1. Таким образом, результаты, полученные здесь на измельчающихся сетках, подтверждают

Е-104

ч

\ / 5х

Рис. 4. Зависимость от времени кинетической энергии возмущений для сетки с шагом Н = 0.025: 1 - а = 0; 2 - а = 0.5; 3 - а = 1; 4 - а = 1.5; 5 - а = 2.

Рис. 5. Зависимость от времени производства кинетической энергии возмущений для сетки с шагом Н = 0.025 (сплошные кривые — расчет по уравнению (19); • — расчет по формуле (20)): 1 - а = 0; 2 - а = 0.5; 3 - а = 1; 4 - а = 1.5; 5 - а = 2.

Таблица 3. Относительные отклонения Де;^(а), %

h а = 0.5 а = 1.0 а = 1.5 а = 2.0

0.1 1.2125 2.5573 5.3571 9.6120

0.05 1.2120 2.5562 5.3548 9.6078

0.025 1.2117 2.5556 5.3536 9.6056

Таблица 4. Относительные отклонения Де)а(Л,), %

h а = 0.0 а = 0.5 а = 1.0 а = 1.5 а = 2.0

0.1 0.06609 0.06690 0.06783 0.06983 0.07312

0.05 0.02203 0.02230 0.02261 0.02328 0.02437

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

Для количественного сравнения вкладов объемной и схемной вязкостей в диссипацию пульсационных величин вычислялись усредненные относительные отклонения

Лел(а)

< E(a, h) > - < E(0, h) >

< E(0, h) >

• 100 %

(21)

при a = 0.5 ... 2,0 и h = const, а также

< E(a, h) > - < E(a; 0.025) >

Ae,a(h)

100%

(22)

< E(a; 0.025) >

для h = 0.1; 0.05 и a = const. Здесь угловые скобки < > означают усреднение по времени

©

< <^(a, h) >= в 1 ^(t, a, h) dt.

(23)

Усреднение проводилось на интервале в = 5. Отметим, что осредненные в соответствии с (23) интегральные характеристики проявляют тенденцию к сходимости по к с порядком 0(к1,5). Результаты расчетов по соотношению (21) сведены в табл. 3.

Как следует из данных табл. 3, при изменении коэффициента объемной вязкости среднее относительное снижение кинетической энергии возмущений приближается к 10 %. При этом изменения Де,^(а) в зависимости от шага сетки к происходят лишь в четвертом знаке после запятой.

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

Заключение

На основе модельной задачи численно исследовано влияние объемной вязкости на нелинейное взаимодействие вихревого возмущения конечной амплитуды с несущим сдвиговым

потоком. Диапазон параметров модельного течения соответствовал реальным значениям для молекулярных газов. Во всех расчетах, выполненных на последовательности вложенных сеток с шагами h = 0.1; 0.05; 0.025; 0.0125, наблюдалась сходимость, в частности, к тестовым точным решениям, с порядком не ниже первого. Полученные результаты позволили количественно выделить вклад схемной вязкости в моделируемые физические дисси-пативные эффекты. Показано, что погрешности, вносимые схемной вязкостью в расчетные пульсационные характеристики, не превышают 0.1 %. Вместе с тем вклад объемной вязкости в диссипацию возмущений составляет до 10 %. Таким образом, подтверждены выводы работы [7], полученные для данной модельной задачи в более широком диапазоне чисел Маха и Рейнольдса, но на более грубой сетке.

В совокупности результаты [7] и настоящей работы дают основание для более детального исследования эффекта объемной вязкости (релаксационного процесса) на основе более совершенных моделей.

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

[1] ЛЕонтович М.А. Замечание к теории поглощения звука в газах jj Журн. экспер. и теорет. физики. 193б. Т. б, вып. б. С. Бб1-Б7б.

[2] A discussion on the firth and second viscosities jj Proc. Roy. Soc. A. 19Б4. Vol. 22б. P. 1-б9.

[3] Кочин Н.Е., КивЕль И.А., Розе Н.В. Теоретическая гидромеханика. Ч. II. М.: Физматгиз, 19б3.

[4] Nerushev A., Novopashin S. Rotational relaxation and transition to turbulence jj Phys. Lett. 1997. Vol. A232. P. 243-24Б.

[Б] Bertolotti F.B. The influence of rotational and vibrational energy relaxation on boundary-layer stability jj J. Fluids Mech. 1998. Vol. 372. P. 93-118.

[6] Григорьев Ю.Н., Ершов И.В. К вопросу о влиянии вращательной релаксации на ламинарно-турбулентный переход jj Тез. докл. Юбилейной науч. конф., посвященной 40-летию Ин-та механики МГУ. Москва, 22-2б ноября 1999. М.: МГУ, 1999. С. бБ-бб.

[7] Григорьев Ю.Н., Ершов И.В. Подавление вихревых возмущений релаксационным процессом в течениях возбужденного молекулярного газа jj Журн. прикл. механики и техн. физики. 2003. № 4. С. 22-34.

[8] Browand F.K., Chin-MinG-Ho. The mixing layer: an example of quasi two-dimensional turbulence jj J. Mecanique Teor. et Appl. 1983. Numero Spec. P. 99-120.

[9] Теория турбулентных струй j Г.Н. Абрамович, Т.А. Гиршович, С.Ю. Крашенинников и др. j Под ред. Г.Н. Абрамовича. М.: Наука, 1984.

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

[11] Корн Г., Корн Е. Справочник по математике для научных работников и инженеров j Под ред. И.Г. Абрамовича. М: Наука, 1970.

Поступила в редакцию 7 февраля 2006 г.