CC BY
2Труды МАИ. Выпуск № 94 http://trudymai.ru/
УДК 533.6.011
Сшивание кинетической и гидродинамической моделей на примере
течения Куэтта
Березко М.Э.*, Никитченко Ю.А.**, Тихоновец А.В.***
Московский авиационный институт (национальный исследовательский университет), МАИ, Волоколамское шоссе, 4, Москва, A-80, ГСП-3, 125993, Россия
*e-mail: jose151294@gmail.com **e-mail: nikitchenko7@yandex.ru ***e-mail: tikhalena@gmail.com
Аннотация
Предложена физико-математическая модель течения газа в пристеночной области. В пределах кнудсеновского слоя используется модельное кинетическое уравнение, в остальной расчетной области - модель Навье-Стокса-Фурье (НСФ). В области сшивания моделей восстанавливается аппроксимирующая функция распределения, представляющая собой разложение локально-равновесной функции Максвелла по степеням тепловой скорости. Параметры разложения (неравновесные напряжения и тепловые потоки) определяются в приближении модели НСФ.
На примере решения тестовой задачи о течении Куэтта показано, что разработанная модель обладает удовлетворительной точностью описания поля течения в широком интервале чисел Кнудсена и Маха.
Ключевые слова: граничные условия, модель Навье-Стокса-Фурье, кинетическое уравнение, сшивание моделей, течение Куэтта.
Введение
Проблема постановки граничных условий на твердой поверхности для уравнений вязкого теплопроводного газа становится весьма актуальной в случае гиперзвуковых и умеренно разреженных течений. Постановка граничных условий на уровне газодинамических переменных приводит к существенным погрешностям при относительно больших значениях чисел Маха (М) и Кнудсена (£п). Физически адекватная формулировка граничных условий на активных (поглощающих или выделяющих газ) поверхностях невозможна без учета молекулярных процессов.
При описании течения на уровне систем моментных уравнений, порядок которых выше порядка системы Навье-Стокса-Фурье (НСФ), помимо проблемы формулировки граничных условий возникает проблема определения достаточного количества таких условий.
Кинетические модели, описывающие течение газа на молекулярном уровне, лишены указанных недостатков. Они позволяют выставлять физически адекватные условия на поверхностях с различными свойствами для любого режима течения. Поскольку граничные условия выставляются для функции распределения молекул по скоростям (весовой функции), значения любого момента этой функции в граничной точке определяются элементарно и единственным образом.
_^
Однако, кинетический расчет поля течения сложной геометрии при Кп < 10 слишком неэкономичен в смысле количества вычислительных операций и необходимого объема памяти компьютера. По этой причине в практических
приложениях для плотных и умеренно плотных газов кинетические модели не используются.
Граничные условия для систем моментных уравнений различного порядка на базе молекулярно-кинетической теории, по-видимому, впервые были сформулированы Грэдом [1]. Дальнейшее развитие эти исследования получили в работах [2, 3, 4].
В граничной точке по параметрам модели, используемой для описания течения, восстанавливалась аппроксимирующая весовая функция падающих на поверхность молекул. Вычисляются потоки массы, импульса, энергии и других моментов на обтекаемую поверхность. Принимается закон взаимодействия молекул с поверхностью. Весовая функция отраженных молекул восстанавливалась с учетом принятого закона и вычисленных потоков.
Аналогичный подход применялся для условий скольжения погранслоя в умеренно разреженных газах [5, 6].
Характерной особенностью известных методов формулирования граничных условий на твердой поверхности является то, что аппроксимирующая весовая функция восстанавливается в граничной точке. В большинстве моделей граничных условий аппроксимирующая функция падающих на поверхность молекул представляет собой разложение локально-максвелловской функции по степеням тепловой скорости.
В условиях течения высокой динамической неравновесности такой вид аппроксимирующей функции не адекватен реальным физическим процессам. Это
связано с тем, что коэффициенты разложения (неравновесные напряжения, тепловые потоки и другие моменты) перестают быть малыми величинами. В периферийных областях пространства скоростей аппроксимирующая весовая функция приобретает отрицательные значения.
Как показывают экспериментальные и теоретические исследования (см., например, [7, 8]) размер высоконеравновесных пристеночных областей ограничен несколькими длинами свободного пробега молекулы. За пределами этих областей указанное разложение весовой функции вполне обосновано.
Цель настоящей работы - разработка физико-математической модели течения в пристеночном слое. Разрабатываемая модель должна описывать течение в пределах кнудсеновского слоя на кинетическом уровне с соответствующими граничными условиями на твердой поверхности. В остальной области течения должна использоваться какая-либо гидродинамическая (моментная) модель, например, модель НСФ (система моментных уравнений неполного второго порядка).
Такая компоновка модели позволит:
• выставить физически адекватные условия на твердой поверхности;
• обеспечить приемлемую для практических реализаций экономичность в широком интервале чисел Маха и Кнудсена.
Последняя особенность модели связана с области, измеряемый в длинах пробега молекулы, слабо зависит от числа Маха.
тем, что размер кинетической не зависит от числа Кнудсена и
1. Гидродинамическая и кинетическая модели
В качестве гидродинамической модели выбрана модель НСФ. В общем случае эта модель может быть представлена следующей системой уравнений [9] (повторяющиеся греческие подстрочные индексы подразумевают суммирование от 1 до 3):
дР | дриа = 0
Ы дха
ди1 ди1 1 дР1а .
- + и„ —- +--— = 0
дг
дТ_
а
+ иг
дха Р дха дТ
. (1)
дхг
+
{г- 1)ТС
ди
ар
а , 1 д^а =0
+
дхр еур дхс
Р-/ =-М
Г ди1 ди/^
дх,- дxi V / 1 у
* 2
+ 8 1 ^ и - 3
с
5 - 3^
1 -
V 2
Н
ди
а .
;дха
(2)
ср дТ
У = — и — 1 Рг дх.
В этих выражениях:
р, и 1, T - плотность, скорость и термодинамическая температура газа;
Р.. = р.. + 8..р - тензор напряжений, Т / =
и и -
рК
удельные напряжения;
р, ргу - давление и тензор неравновесных напряжений;
- вектор теплового потока;
ии = ) - коэффициент вязкости, ^ - свободный параметр модели;
И - свободный параметр модели, имеющий смысл отношения времени поступательно-вращательной релаксации к времени поступательно-поступательной релаксации энергии молекул, в настоящей работе принять к = 5;
с , ср у, Я - изохорная и изобарная теплоемкости, показатель адиабаты,
удельная газовая постоянная.
Граничные условия на твердой поверхности модели НСФ формулируются как условия скольжения погранслоя [2]. В разрабатываемой модели течения такие граничные условия не требуются.
В качестве кинетической модели рассматривается модельное кинетическое уравнение многоатомных газов, аналогично работе [10]:
^=(3)
д1 а дха Тр
В этом выражении:
^ - молекулярная скорость (координата пространства скоростей);
и
т р =--время релаксации напряжений;
р Р
/ = /х1 ,х2,- весовая функция, фазовое пространство, которой
дополнено подпространством вращательной энергии молекулы е.
После формального интегрирования по пространству вращательных энергий (3) может быть приведено к двум уравнениям:
д_ дг
/ /г
д
дх,
/г
/г
т,
/+- /г
/г /г
Здесь:
/г = I !г = [£/А£;
/+
п
{2кят;}
'/2
ехр
саса
2ЯТ+
1 +
РаСс
3р{кТ+ + )* V 5КТг
саса _ ^
/++ = ^^ кТ+ /++;
2{Г-1)
Т+ += Т +
5 - 3^
2
1 - 1)(Тг - Тг); Тг+ = Т - 3{Г -1/1 - - Тг);
сг = ^ - и - тепловая скорость движения молекул; величина 1/7 представляет долю неупругих столкновений по отношению к общему числу столкновений молекул.
Макроскопические параметры газа определяются следующим образом:
п = I/А; Р =т01^/¡А; Т = Раа/3РК; п^Тг = 1;
Гт0С Г!
Рг
2
здесь: т 0 и k - масса молекулы и постоянная Больцмана. Для формулировки граничных условий на твердой поверхности достаточно задать закон взаимодействия молекул газа с поверхностью. В настоящей работе
использован закон диффузного отражения с полной аккомодацией энергии. Этот
1
закон вполне адекватен реальным физическим процессам, протекающим на шероховатых поверхностях при умеренных энергиях (< 1эв) движения молекул [3].
2. Метод сшивания моделей
Разработанная модель течения (модель КИН_НСФ) предполагает разбиение поля течения на две области - кинетическую и гидродинамическую. Кинетическая область непосредственно граничит с обтекаемой поверхностью. Ее протяженность в направлении нормали к поверхности составляет несколько длин свободного пробега молекул X (кнудсеновский слой). Остальное поле течения занято гидродинамической областью.
На смежных границах указанных областей расположена область сшивания кинетической и гидродинамической моделей. В кинетической и гидродинамической областях соответствующие модели течения реализуются традиционным для них образом. В области сшивания численное решение моделей специфично.
Известно [11], что кинетические модели более информативны, чем гидродинамические. Передача информации в виде значений р, щ, Т, р от кинетической модели к гидродинамической не вызывает затруднений. Указанные газодинамические параметры являются моментами весовой функции, определяемой кинетической моделью.
Обратный информационный обмен неоднозначен. По газодинамическим параметрам, определенным гидродинамической моделью, принципиально невозможно восстановить весовую функцию в динамически неравновесном течении.
В условиях, близких к равновесным, возможны адекватные аппроксимации функции. Отметим, что модель НСФ теоретически обоснована только в слабо неравновесных течениях (модель первого приближения процедуры Чепмена-Энскога). Этой модели соответствует весовая функция, часто называемая "навье-стоксовой" [2].
Если область сшивания моделей находится в слабонеравновесной области течения, то передача информации от модели НСФ к кинетической модели посредством весовой функции типа навье-стоксовой будет вполне обоснована. Модель течения и аппроксимирующая весовая функция будут соответствовать одному порядку приближения. Пригодность модели течения для описания процессов в рассматриваемой области будет означать пригодность аппроксимирующей функции для передачи информации в кинетическую модель.
Известно [2], что профиль скорости в погранслое меняет логарифмическую форму на линейную на расстоянии от поверхности. Линейный профиль
скорости характерен для модели НСФ. На этом основании сделана приведенная выше оценка толщины кинетической области.
В настоящей работе принята аппроксимирующая весовая функция /А, аналогичная навье-стоксовой. В случае многоатомных газов интегралы этой функции по пространству вращательных энергий приобретают вид:
Г 2 \ У\
/Лг = /Мг
1 +
ртят
1 т
~^рарсаср +
С 1
5ЯТ
<РаСа
V У ±"г V Ч^^ у ])
/Аг - кТг
5 - Ъу 2(у -1)
fAt + /ш
®аса
ртКТг
; /ш1 -
п
г У
(¿жЯТ У2
ехр
( с2 Л
2 КТ,
Параметры функций /А(, /А1 определяются в приближении модели НСФ [11]:
Т - Т - ^У ; Тг - Т + (у - 1) ЬК^а ; рт - рЩ ;
Ъ Кр дх
а
Кр дх
а
Р[
Л
Гди- дм,- 2 ди Л +__^ — ^
дх, дх,- у Ъ дх
У
15 дТ
К -—~тКл—; ^ - — 4 дх,
/
су - - К У 2
л
дТ
У дх1
л
На рис.1 показана схема расчетной области в предположении, что для численного решения задачи используется конечно-разностный метод. Область сшивания моделей обозначена символом С с тремя перечеркнутыми узлами.
Рис.1. Схема расчетной области. А - кинетическая область, В
гидродинамическая область, С - область сшивания.
В узлах кинетической области (область А, темные точки) решается
кинетическое модельное уравнение (4). Параметры р, щ, Т, р, определенные в
верхней точке области А, служат краевыми значениями для модели НСФ. Эта
10
модель решается в гидродинамической области (область В, светлые точки). Аппроксимирующие весовые функции восстанавливаются в двух верхних узлах области сшивания.
Количество таких узлов зависит от выбора конечно-разностного шаблона кинетической модели. На рис. 1 предполагается аппроксимация производной по трем точкам, что требует двух краевых узлов. Для модели НСФ достаточно одного краевого узла, если для решения поставлена краевая задача, реализуемая на трехдиагональной ленточной матрице.
На обтекаемой поверхности, в нижнем узле области А, выставляются кинетические граничные условия в их традиционной форме.
3. Тестовая задача
В качестве тестовой задачи выбрано плоское течение Куэтта. Это одномерное течение хорошо изучено как экспериментально, так и теоретически [12].
Схема течения и используемая система координат показаны на рис.2. Две плоские бесконечные поверхности движутся в противоположных направлениях оси OX с относительной скоростью uW. Температура поверхностей - TW. До начала движения плотность газа и концентрация молекул имели значения р0 и п0.
Т 1 1 \Л/ 1 1У 0.5 -
X
0
Т\л/ 1 — 0.5 Цд,
Рис. 2. Схема тестовой задачи.
Средняя длина свободного пробега молекулы в невозмущенном (неподвижном) газе:
^ = Мо
Числа Маха и Кнудсена определены как
М = ; Кп = Мо 1,
здесь Мо - коэффициент вязкости, рассчитанный по температуре Тш . Варьирование числами М и Кп позволяет изучать течения различной степени неравновесности. Величина касательного напряжения рху может приближаться к
величине давления р.
Особенностью течения Куэтта является обращение уравнения сохранения массы (первое уравнение системы (1)) в тривиальное тождество. Поле плотности определяется по уравнению сохранения импульса, записанному для иу. В
стационарном течении иу =0 и закон сохранения импульса сводится к
соотношению
dP
-PZ- = 0. -y
В рассматриваемом течении в приближении модели НСФ pyy = 0 (см. (2)), следовательно Pyy = p = pRT = const. Очевидное соотношение
2 L 2
- Jpdy = po (5)
L o
может быть преобразовано следующим образом:
-Lfpdy = - J PRTdy = -pRT?±_
L J L J RT L J RT
Отсюда следует:
p(y)= LP . • (6)
2T (y )J
dy
Т (у )
Поле плотности определяется локально алгебраическим уравнением (6). В отличие от дифференциальных уравнений, определяющих их и Т, (6) не связано явно с фиксированными параметрами uw и Tw на границах вычислительной области. Это приводит к накоплению вычислительной погрешности в процессе итерационного решения (1). Аналогичная проблема возникает при решении (4).
Для коррекции значений плотности на каждой итерации используется соотношение (5). Корректированное значение плотности (р') определяется как
Р'-^р- (7)
ь о 13
Такая коррекция носит приближенный характер, так как требует постоянства только средней по полю плотности. При надлежащем выборе схемы решения вычислительная погрешность в пределах одной итерации имеет место только в последних знаках мантиссы числа. При проведении коррекции (7) на каждой итерации неточности коррекции практически не влияют на установившиеся значения параметров течения, определенных кинетической и НСФ моделями.
В модели КИН_НСФ влияние коррекции на установившееся решение несколько больше. Это связано с тем, что кинетическая модель имеет большую вычислительную погрешность, чем модель НСФ. Дополнительную погрешность вносит численное интегрирование по пространству скоростей. В результате коррекция (7) оказывает завышенное влияние на гидродинамическую область и заниженное - на кинетическую.
Отметим, что коррекция плотности вызвана спецификой течения Куэтта. В других задачах, например в задачах внешнего обтекания, фиксированные значения плотности (рж) присутствует на границе вычислительной области. Проблема накопления вычислительной погрешности не возникает.
4. Результаты тестовых расчетов
На начальном этапе тестирования моделей выбирались свободные параметры я, И и Ъ из условия наилучшего соответствия расчетных и экспериментальных
коэффициентов трения в криптоне при М=0.88 и в воздухе при М=1.4 [12]. Интервал
_^
чисел Кнудсена - Кп = 10 ...1. Для модели НСФ принято: я = 0.7, к = 5. Для
кинетической модели: ^ = 1, Z = 5. Эти же значения использовались в гидродинамической и кинетической областях модели КИН_НСФ. Область сшивания этой модели находилась на расстоянии 2Х от твердой поверхности.
Результаты расчетов показали, что наиболее наглядно расхождения профилей наблюдается на профилях температуры. На рис. 3 показаны профили температуры двухатомного газа при М = 1 и Кп = 0.01. Здесь и на следующих рисунках профили кинетической модели изображены сплошной линией, модели НСФ - мелким пунктиром, модели КИН_НСФ - крупным пунктиром.
У
Рис. 3. Профили температуры, М = 1, Кп = 0.01. Сплошная линия -кинетическая модель; мелкий пунктир - модель НСФ; крупный пунктир - модель КИН НСФ.
При данном числе Кнудсена модель НСФ с граничными условиями скольжения погранслоя теоретически обоснована. Профили всех моделей практически совпадают. Наибольшее отклонение от кинетического имеет профиль
модели КИН_НСФ. Это отклонение связано с указанными выше погрешностями коррекции концентрации. Рис.3 позволяет оценить влияние этих погрешностей.
При данных параметрах течения модель КИН_НСФ потребляла в несколько десятков раз меньше процессорного времени и оперативной памяти, чем кинетическая модель.
На рис.4 представлены результаты расчета при М = 1 и Кп = 0.1. Это значение Кп соответствует переходной области течения, в которой модель НСФ не обоснована теоретически, что и подтверждает рис.4. Модель КИН_НСФ, использующая НСФ в большей части расчетной области, дает существенно лучший профиль, даже с учетом погрешности коррекции плотности.
У
1.01 1.02 1.03 1.04 1.05
Рис. 4. Профили температуры, М = 1, Кп = 0.1.
Графики рис.5 построены при М = 4, Кп = 0.01. Из сравнения с рис.3 видно,
что увеличение числа Маха до гиперзвуковых значений не существенно влияет на
погрешность модели НСФ пока число Кнудсена находится в континуальной области
течения. Отклонение профиля модели КИН_НСФ не сильно отличается от рис.3,
16
если рассматривать это отклонение в относительном измерении.
1 1.2 1.4 1.6 1.8
Рис. 5. Профили температуры, М = 4, Кп = 0.01.
На рис.6 приведены графики для М = 4, Кп = 0.1. Сравнение с рис.4 показывает существенные отклонение профиля модели НСФ в верхней части рисунка, соответствующей пристеночной области. Это свидетельствует о том, что в переходной области течения при больших скоростях граничные условия скольжения погранслоя престают адекватно описывать процесс взаимодействия газа с твердой поверхность.
у
1.2 1.4 1.6 1.8
Рис. 6. Профили температуры, М = 4, Кп = 0.1. На рис.7 показана зависимость среднеквадратического отклонения (а) профилей температуры (Т/Тф) моделей НСФ и КИН_НСФ от профиля кинетической модели в зависимости от чисел М и Кп.
Рис. 7. График среднеквадратичного отклонения профилей температуры в зависимости от чисел Маха и Кнудсена.
Характер зависимостей вполне ожидаем. Увеличение чисел М и Кп означает повышение степени неравновесности течения, оцениваемой как произведение МКп [11]. Модель НСФ теоретически обоснована только в слабонеравновесных течениях, что и подтверждают кривые на рис.7. Отклонения модели КИН_НСФ существенно меньше.
Обращает на себя внимание кривая модели КИН_НСФ при значениях
Кп = 0.01, отвечающих минимальной степени неравновесности в рассмотренном
интервале чисел Кп. Эта кривая близка к соответствующей кривой модели НСФ, и даже несколько превосходит ее. К тому же, что эта кривая проходит существенно выше кривых модели КИН_НСФ, соответствующих большей степени неравновесности.
По-видимому, такое поведение модели КИН_НСФ объясняется неточностью коррекции концентрации, так как при этом числе Кп обе сшитые модели достаточно хорошо описывают течение. Уместно напомнить, что коррекция концентрации является специфической процедурой решения задачи о течении Куэтта.
Заключение
Результаты расчетов показали, что модель КИН_НСФ ненамного уступает в точности кинетической модели и существенно превосходит ее в экономичности. При Кп = 0.01 модель КИН_НСФ требовала в несколько десятков раз меньше процессорного времени, чем кинетическая модель.
При описании течений плотных газов экономичность модели КИН_НСФ не зависит от числа Кп и слабо зависит от числа М. Таким образом, разработанная модель течения в отношении адекватности граничных условий наследует свойства кинетической модели, а в отношении экономичности - свойства модели НСФ.
Библиографический список
1. Грэд Г. О кинетической теории разреженных газов: перевод с английского // Механика. 1952. №4. С. 71-97.
2. Коган М.Н. Динамика разреженного газа. - М.: Наука, 1967. - 440 с.
3. Баранцев Р.Г. Взаимодействие разреженных газов с обтекаемыми поверхностями. - М.: Наука, 1975. - 343 с.
4. Сакабеков A.C. Начально-краевые задачи для систем моментных уравнений Больцмана в произвольном приближении // Математический сборник. 1992. Т. 183. № 9. С. 67-88.
5. Латышев А.В., Юшканов А.А. Моментные граничные условия в задачах скольжения разреженного газа // Известия РАН. Механика жидкости и газа. 2004. № 2. С. 193-208.
6. Никитченко Ю.А. Феноменологическая модель граничных условий на твердой поверхности // Вестник Московского авиационного института. 2012. Т. 19. №3. С. 514.
7. Becker M., Boyland D.E. Flow field and surface pressure measurements in the fully merged and transition flow regimes on a cooled sharp flat plate // Rarefied Gas Dynamics, Suppl. 4. 1967. Vol. 2, pp. 993-1014.
8. Tannehill J.C. ,Mohling R.A., Rakich J.V. Numerical computation of the hypersonic rarefied flow near the sharp leading edge of a flat plate // AIAA Paper. 1973. №73-200, pp. 1-13.
9. Никитченко Ю. А. Модели первого приближения для неравновесных течений многоатомных газов // Труды МАИ. 2014. №77. URL: http://www.mai.ru/science/trudy/published.php?ID=52938
10. Glinkina V.S., Nikitchenko Yu.A., Popov S.A., Ryzhov Yu.A. Drag Coefficient of an Absorbing Plate Set Transverse to a Flow // Fluid Dynamics. 2016. Vol. 51. No. 6. pp. 791 - 798.
11. Никитченко Ю.А. Модели неравновесных течений - М.: Изд-во МАИ, 2013. -160 с.
12. Кошмаров Ю.А., Рыжов Ю.А. Прикладная динамика разреженного газа: - М.: Машиностроение, 1977. - 184 с.
