CC BY
130
УДК 519.6
Ю. И. Димитриенко, А. А. Захаров, А.С. Абакумов, М. Н. Коряков, Е. К. С ы з д ы к о в
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ГАЗОВЫХ ПОТОКОВ В КАНАЛАХ ВОЗДУХОЗАБОРНИКОВ НА ОСНОВЕ УРАВНЕНИЙ НАВЬЕ-СТОКСА
Предложен алгоритм численного решения уравнений Навье-Стокса в адаптивных координатах для областей сложной формы, основанный на методе ленточно-адаптивных сеток и схеме TVD 2-го порядка точности. Представлены результаты численного тестирования алгоритма задачи о затекании сверхзвукового потока в канал со ступенькой и течении в канале сверхзвукового воздухозаборника. Проведено сравнение режимов течения для идеального и вязкого газов, которое показало, что для вязкого газа значения параметров потока — давления, плотности и температуры на стенках — существенно превышают соответствующие значения для идеального газа.
E-mail: dimit.bmstu@gmail.com
Ключевые слова: газовая динамика, уравнения Навье-Стокса, численные
методы, ленточно-адаптивные сетки, схемы TVD.
В настоящее время для численного моделирования процессов газовой динамики применяют как коммерческие расчетные пакеты, так и собственные разработки компаний, занимающихся проектированием изделий, разработки университетов и научно-исследовательских институтов. По-видимому, основная причина, почему коммерческие пакеты не вытеснили собственные разработки компаний и университетов, заключается в том, что универсальные методы вычислений, заложенные в коммерческие продукты, не гарантируют точности вычислений при решении сложных задач, которые предварительно не тестировались их разработчиками. Хорошо известно, что ни один из современных вычислительных методов не обладает абсолютными преимуществами по качеству получаемого решения и не является универсальным, пригодным для широкого набора задач газовой динамики.
Узкоспециализированные программы часто позволяют добиваться превосходства над более универсальными пакетами типа ANSYS FLUENT, STAR CD, ANSYS CFX и т.п. как по быстродействию, точности решения и минимизации требуемых ресурсов компьютера, так и по удобству использования программ и их настройки.
Численное решение точных уравнений Навье-Стокса для областей сложной геометрической формы относится к числу задач, достаточно трудных для коммерческих программных пакетов (несмотря на формально заявляемые ими возможности решения "любых" задач). В настоящей статье представлены результаты работ по развитию численного метода ленточно-адаптивных сеток [1] для системы уравнений
Навье-Стокса в областях сложной формы типа каналов сверхзвуковых воздухозаборников прямоточных воздушно-реактивных двигателей. Эти работы продолжают ранее выполненные исследования [1] по разработке метода ленточно-адаптивных сеток для уравнений Эйлера.
Система уравнений Навье-Стокса для вязкого теплопроводного газа содержит уравнение неразрывности, три уравнения движения и уравнение энергии в бескоординатной форме имеет следующий вид: др „
ж+= 0;
Яру
' + У-(ру <8> у + рЕ - Т*) = 0; (1)
dt
dpE
+ V ■ ((pE + p) v - Tv ■ v + q) = 0.
дг
К этим уравнениям присоединяются определяющие соотношения вязкого теплопроводного газа
р = Ярв; е = с„в; Е = е + |у|2 /2; (2)
Т = р1(У-у)Е + + У®ут); (3)
я = -лvе, (4)
где Т — тензор вязких напряжений; я — вектор потока теплоты; — коэффициенты вязкости газа; Л — коэффициент теплопроводности газа; р — плотность; Е — плотность полной энергии газа; е — плотность внутренней энергии; су — теплоемкость при постоянном объеме; в -температура; у — вектор скорости; |у |2 = у ■ у — квадрат модуля скорости; р — давление; Я — газовая постоянная; Е — метрический тензор; V - набла-оператор.
Уравнения Навье-Стокса в адаптивных координатах. Рассмотрим три типа координат: декартовы X, адаптивные криволинейные координаты X^, которые согласованы с границей рассматриваемой геометрической области течения потока, и ортогональные криволинейные координаты X,3 (например, цилиндрические). Тогда система уравнений Навье-Стокса (1) может быть записана в символическом виде в адаптивных координатах X^ [1] в недивергентном виде:
д ^р Л- ( а ^ =0
дг ¿1Ра дх^ На ) 0'
d3 ^ +
+£ НВ?* - pi+HE (^ - P1 p*)
/-V /4—1 \ \ / /
+ E p УНГ^ + p) ) =
31 ^ адХ1 у на V р
= т (+V м«^ (р>
" ^ Ра дХ5 нЛ на дХ- + ^ м нддХ- + 1
Здесь обозначено: Vа — компоненты вектора скорости в физическом базисе г ъ координат X°; д'а/3 = ^ — метрическая матрица
этого базиса;
Л? - НИИ • P - 1 dHa X 1 1Нв X . - Н1Н2Нз. Гв„ - дщ 9X0 - И ÖXY•
1Xl 1xl p'
pj - ^ • Q'j - o^j • на -yfgaa • - •
(6)
p _ _ _ 1 /р я R 3 p
rpaß _ jMj-aßklRlk _ \ л j^aßsu 1 I pS __^ , \ л гш
- - Hf\ £ ÖX'S ^
J£
£,U>=1 ь ^ 7=1
kl I ,, t xik <
Mljkl - ^18ljSkl + p,2(ölkSjl + 5ÜSjk). Запишем систему (1) в полных дифференциалах [1]:
1 i + f р-А. 1 ^ - у P.JL (^)- (7)
PV9\dtp + ^P'ÖXIР 2-,p«dx\ На г (7)
=1 =1
nt'dvY ^ 1 (PS (р ÖVY r Öp \ р р -
IT + ^ WApa[V*1X + ^!Xs) + v«vsГ"
a,s=1
¿( pa м На+На р^>
/—(de ^ - / va de д f л/gv'-
d + £ Pa H dx+
dt ^ Ha dXs "дЧ Ha
a,s=1 4
Л ps d f^g \P*a de \ r- „ adXS dX^ + V ^ ,
где
ЕМда/Зеш / дЬо - \ дV - -
_- -в + V Гв Р*+ V Гш
НаН£ Уа дХ* + °5 г / V дХа + ° "г "
— функция диссипации.
Граничные условия к системе уравнений (7) в адаптивных координатах имеют вид:
дв
V'1 = 0, -ЛдХХ- П = д£ (8)
— на жесткой стенке;
р = ре, У' = VI, в = ве (9)
— на границе входа потока (здесь Vге, ве — заданные значения скорости и температуры);
• д 7игп дв .
П — = 0, — П = 0 (10)
дХ? дХг 4 у
— на сверхзвуковой границе выхода потока;
др
дХ
-иг = 0; vг n г = 0;
(11)
dvmт дв
п __lam = о a = 12- ——П = 07 д> дв i n n^^T- = 0, ——- П = 0 дХ- дХг
— на поверхности симметрии. Здесь ni — компоненты вектора нормали и тam — компоненты векторов в касательной плоскости. Начальные условия в адаптивных координатах имеют вид
t = 0 : р (0,Хг) = р0, v4 (0,Х-) =0, E (0,Хг) = cvв0. (12)
Использование метода ленточных адаптивных сеток. Генерация геометрически-адаптивных сеток основана на введении специальной криволинейной системы координат Х- (жг), согласованной с границей рассматриваемой области. Сетка имеет глобальную сквозную нумерацию узлов и называется ленточной [1]. В трехмерном случае для
каждого узла можно алгоритмически найти его 6 соседей, за исключением граничных и угловых точек.
Система уравнений (5) имеет второй порядок производных и для ее решения был реализован модифицированный метод ЛАС. Модификация осуществляется на основе метода расщепления по физическим процессам, суть которого в следующем: на каждом временном шаге вычисления разбиваются на два этапа.
На первом этапе решается система в виде (5) без вязкости и теплопроводности:
д^Р , ^ А д (
dt +2. ра dj ж pv°> =
а= 1
dt( +
d
{Vg'Pz
+ Е (ра dj H. +наRß =o <13>
dt
^ d (VF-аЛ г , V^
^+± й и € р-( Е+Р)) =0.
Для ее численного решения применяется метод ЛАС на основе схемы ТУБ [2, 3], значения же в граничных узлах вычисляются с помощью разностной аппроксимации граничных условий так же, как и для идеального газа.
На втором этапе решается параболическая часть системы в виде (7) без конвективных членов и без уравнения неразрывности:
Р^ £-¿(* д^ ( £ ^ ) + £ ГТ!,);
(14)
a,s=1
Г-,de Л pS d /V97Ара дв \ + г- „ PV9'dt = ра^-WaHöX«) •
Для ее решения применяется метод, основанный на расщеплении дифференциальных операторов по координатным направлениям. В этих целях система (14) записывается в векторном виде
ргдЦ = Лиз + (15)
где дифференциальный оператор имеет вид
aus - pk
д
5 P Втт
EAlj pП ÖU Z
4 SZpj ÖXn
Z=1
+ Bjz pj QXr
ÖP 5 ^ U Z + У" Csz Uz, (16)
Z=1
Z=1
а координатные столбцы Us и Fs определяются так:
Us -
(р л (0\
P 2 • F s - 0
P 3 5 5 0
\pv в)
Ненулевые матрицы Aj имеют следующий вид:
Aj+ilß+i - Vg'M^, а, ß -1, 2,3:
Aj+i -\fg' MMmj «m, ß -1, 2, 3i
(17)
A5j5 - \fg' Xg'lj, Mjk - Mmjspp'p'köPl.
Для численного решения полученной системы квазилинейных параболических уравнений (15), (16) применим метод расщепления по координатным направлениям, используя для этого экономичные разностные схемы, пригодные для уравнений с переменными коэффи-
диз . ^
циентами. Аппроксимируя операторы и лиз соответствующими разностными выражениями, в результате получаем четырехшаговую
разностную схему:
Um+1/4 - Й
At Л т+1/4 - YAlUs ,
m I um+1/2 — um+1/4
U'm+3/4 _ u"m+1/2
pm (и'+1 - и'+3/4) - A (Ai + Л2 + Лз)
+(Aq Usm + Fm )At,
- AA2Usm+1/2,
2
At Л TTm+3/4
- TA3Us ,
At
T
(18)
3) U'+3/4 +
где рт, ит, — значения функций после первого этапа метода расщепления по физическим процессам, полученные с помощью решения
" гл /1\ т7-т+1/4 7-7-т+1/2 пт+3/4 ттгп+1
системы уравнений (14), а и , и , и , ит+1 — значения функций на промежуточных шагах разностной схемы. Линейные дифференциальные операторы Л1 соответствуют координатному на-
m
m
правлению вдоль адаптивнои координаты
Л ? — ра(V A j Ра+ V B Ра•
ЛаТ s Pi QXa dXa / ^ ^szdXa '
z=1 ' z=1
ЛТ - E ра ¿ (¿ ASZpe Ш + ¿ Csz ?. <19)
а,в=1 z=1 ' z=1
а=в
Первые три разностных уравнения в системе <18) решаются метода-
т т-ш+1/4 j-T-m+1/2 r,-m+3/4
ми прогонки относительно ís , ís , ís соответственно, а последнее уравнение — по явноИ схеме относительно í'm+1-
Задача о затекании сверхзвукового потока в канал со ступенькой. Для исследования влияния вязкости газа на результаты решения была рассмотрена одна из классических задач для численных методов — об обтекании канала со ступенькой сверхзвуковым потоком. Были проведены расчеты для числа Маха в невозмущенном потоке M — 3 (рис. 1, 3-я полоса обложки). Левая граница в рассматриваемой области является границей входа потока, правая — границей выхода потока, верхняя и нижняя границы представляют собой жесткую стенку. Коэффициент Пуассона газа y — 7/5, соотношение высоты и длины канала 1 : 3, безразмерная высота ступеньки, отнесенная к длине и равная 0,2, располагается на расстоянии 0,6 от левой грани. Начальные условия в безразмерном виде: p — 0,05, р — 0,5, v1 — v2 — 0; граничные условия на границе входа потока: p — 1/y, р — 1, v1 — 3, v2 — 0. Задача решалась в декартовой системе координат. Расчеты проводились на суперкомпьютере СКИФ-МГУ "Чебышев".
На рис. 1-4 (3-я полоса обложки) приведены сравнительные результаты расчетов для параметров потока в случае идеального и вязкого газов. Картина течения вязкого газа существенно отличается от соответствующей картины для идеального газа: система скачков, образующаяся в вязком газе, изменяет свое расположение, изгибаясь по потоку, кроме того появляются дополнительные скачки, идущие от начальных точек нижней и верхней стенок канала. Эти скачки обусловлены скачком граничных условий в угловых точках, который неизбежно возникает при точном удовлетворении граничным условиям в уравнениях Навье-Стокса (на твердых стенках скорость нулевая, а на границе входа — максимальная). В окрестности твердой обтекаемой поверхности канала образуется пограничный слой, в котором продольная скорость потока v1 резко меняет значения от нулевых до максимальных. Тестирование численных методов расчета для идеального газа описано в работе [2].
Были проведены исследования влияния числа разбиений сетки на результаты расчета: при крупной сетке 50 х 150 <6501 узел) получалось
Иллюстрации к статье Ю.И. Димитриенко, А.А. Захарова, А.С. Абакумова, М.Н. Корякова, Е.К. Сыздыкова «Численное моделирование газовых потоков в каналах воздухозаборников на основе уравнений Навье-Стокса»
Рис. 1. Распределение набегающей скорости, м/с, в канале со ступенькой для идеального (а) и вязкого (б) газов
(а)
(б)
Рис. 2. Распределение плотности, кг/м3, в канале со ступенькой для идеального (а) и вязкого (б) газов
(а)
(б)
Рис. 3. Распределение давления, Па, в канале со ступенькой для идеального (а) и вязкого (б) газов
(а) (б)
Рис. 4. Распределение температуры, К, в канале со ступенькой для идеального (а) и вязкого (б) газов
V
У
нефизическое решение, при котором происходило запирание потока на входе в канал. Для более мелкой сетки 120x150 (15391 узел) удавалось получить адекватное с физической точки зрения решение, на котором были видны все основные скачки, образующиеся в канале. Дальнейшее измельчение сетки до 200x300 (50901 узел) не приводило к качественному изменению решения, а лишь уменьшало зону размытия скачков и толщину пограничного слоя.
Проведены расчеты для неравномерных сеток, сгущающихся к твердым поверхностям по степенному закону. Было установлено, что использование неравномерных сеток ухудшает качество решения, приводит к исчезновению скачков: решение вблизи твердых поверхностей "забивает" решение в основном канале.
Также исследовалось влияние задания начальных и граничных значений в угловой начальной точке с использованием следующего алгоритма. Вначале проводился расчет течения вязкого потока, обтекающего плоскую пластину. На установившемся режиме этого течения выделялась конечная зона выхода потока, на которой формировался профиль продольной скорости , близкий к параболическому. Эти значения скорости ^1(ж2), плотности р(ж2) и давления р(ж2) выбирались в качестве граничного условия в угловой точке в основной задаче. Они обеспечивали численную гладкость решения, поскольку удовлетворяли условию ^1(0) = р2(0) = 0 и ^1(ж2) = ре, где Ре — значение продольной скорости набегающего потока. Сравнительные результаты расчетов для этого алгоритма и граничных условий с разрывом в угловой точке показаны на рис. 5 (4-я полоса обложки). Общая картина течения при устранении разрыва в угловой точке сохраняется, но происходит некоторое уменьшение угла наклона всех скачков в канале.
Расчет сверхзвукового воздухозаборника проведен для осесим-метричного потока в канале со следующими параметрами на входе и начальными данными:
р^ = 0,1946 кг/м3, г^ = 900 м/с, = 0, р^ = 12084,5 Па;
р0 = 0,1946 кг/м3, рг0 = 0, = 90 м/с, р0 = 12084,5 Па.
На рис. 6, 7 (4-я полоса обложки) приведены результаты расчетов параметров потока в канале и вне воздухозаборника для случая идеального и вязкого газов. Число Маха в невозмущенном набегающем потоке М = 3. При численном решении использована адаптивная система координат X0, а в качестве координат X0 использована цилиндрическая система координат. Число узлов сетки — 41258. При расчете течения идеального газа и на первом этапе расчета вязкого газа (невязкое течение) применялись схемы ТУБ. Результаты сравнения показывают значительное различие в решениях для идеального
Иллюстрации к статье Ю.И. Димитриенко, А.А. Захарова, А.С. Абакумова, М.Н. Корякова, Е.К. Сыздыкова «Численное моделирование газовых потоков в каналах воздухозаборников на основе уравнений Навье-Стокса»
(а) (б)
Рис. 5. Распределение числа Маха в канале со ступенькой для вязкого газа в случае граничных условий с разрывом в угловой точке (а) и без разрыва (б)
(а) (б)
Рис. 6. Распределение плотности в канале воздухозаборника для идеального (а) и вязкого (б) газов
(а) (б)
Рис. 7. Распределение температуры в канале воздухозаборника для идеального (а) и вязкого (б) газов
и вязкого газов: для вязкого газа происходит торможение потока на стенках канала воздухозаборника, что приводит к повышению значений плотности в самом канале (рис. 6, б) и температуры и давления на стенках канала (рис. 7, б). Максимальные значения температуры на стенках для вязкого газа достигают значений 850 K на входном конусе и 919 K на внутренней поверхности обечайки воздухозаборника в окрестности первого косого скачка, в то время как в идеальном газе максимальные значения температуры на стенках составляют 660 K (стенки предполагались теплоизолированными).
Выводы. Предложен алгоритм численного решения задач динамики вязкого теплопроводного газа в адаптивных криволинейных координатах на основе схем TVD 2-го порядка и метода ленточных адаптивных сеток. Проведенное численное для задач о течении сверхзвукового потока в канале со ступенькой и о течении в сверхзвуковом воздухозаборнике показало существенное различие в характере потоков для вязкого и идеального газов: в вязком газе за счет торможения потока в окрестности твердых стенок возрастают значения плотности, давления и температуры газа, изменяются углы наклона скачков уплотнения, а также могут возникать дополнительные скачки, вызванные наличием дополнительных разрывов потока в угловых точках при натекании и резком торможении вязкого газа на твердой поверхности.
Работа выполнена при поддержке гранта Президента Российской Федерации № МК-2498.2011.8.
СПИСОК ЛИТЕРАТУРЫ
1. Димитриенко Ю. И., З а х а р о в А. А. Метод ленточных адаптивных сеток в газовой динамике. - М.: Изд-во НТЦ "Университетский", 2008. - 180 с.
2. Д и м и т р и е н к о Ю. И., Коряков М. Н., Захаров А. А., С ы з д ы к о в Е. К. Развитие метода ленточно-адаптивных сеток на основе схем TVD для решения задач газовой динамики // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. - 2011. - № 2. - С. 87-97.
3. ГильмановА. Н. Методы адаптивных сеток в задачах газовой динамики. - М.: Физматлит, 2000. - 248 с.
Статья поступила в редакцию 7.07.2011
Юрий Иванович Димитриенко родился в 1962 г., окончил МГУ им. М.В. Ломоносова в 1984 г.. Д-р физ.-мат. наук, профессор, заведующий кафедрой "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана, действительный член Академии инженерных наук. Автор более 160 научных работ в области вычислительной механики, нелинейного тензорного анализа, термомеханики композитов, математического моделирования в материаловедении.
Yu.I. Dimitrienko (b.1962) graduated from the Lomonosov Moscow State University in 1984. D. Sc. (Phys.-Math.), professor, head of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Full member of the Russian Academy of Engineering Sciences. Author of more than 160 publications in the field of computational mechanics, nonlinear tensor analysis, thermomechanics of composite materials, mathematical simulation in science of materials.
Михаил Николаевич Коряков родился в 1987 г., в 2010 г. окончил МГТУ им. Н.Э. Баумана. Аспирант кафедры "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана. Автор ряда работ в области вычислительной газодинамики.
M.N. Koryakov (b. 1987) graduated from the Bauman Moscow State Technical University in 2010. Post-graduate of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Author of several publications in the field of computational gas dynamics.
Андрей Алексеевич Захаров родился в 1982 г., в 2005 г. окончил МГТУ им. Н.Э. Баумана. Канд. физ.-мат. наук, доцент кафедры "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана. Автор 20 научных работ в области вычислительной газодинамики.
A.A. Zakharov (b. 1982) graduated from the Bauman Moscow State Technical University in 2005. Ph. D. (Phys.-Math.), assoc. professor of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Author of 20 publications in the field of computational gas dynamics.
Александр Сергеевич Абакумов родился в 1987 г., окончил МГТУ им. Н.Э. Баумана в 2010 г.. Инженер НИЧ НУК ФН МГТУ им. Н.Э. Баумана. Автор ряда работ в области вычислительной газодинамики.
A.S. Abakumov (b. 1987) graduated from the Bauman Moscow State Technical University in 2010. Engineer of Research Division of "Fundamental Sciences" scientific-educational complex of the Bauman Moscow State Technical University. Author of a number of publications in the field of computational fluid dynamics.
Елтуган Кимашевич Сыздыков родился в 1956 г., в 1979 г. окончил КГТУ им. А.Н. Туполева. Канд. техн. наук, заместитель генерального директора ОАО ГосМКБ "Радуга" им. А.Я. Березняка. Автор более 50 научных работ в области механики аэрокосмических конструкций.
Ye.K. Syzdykov (b. 1956) graduated from the Kazan State Technical University n.a. A.N. Tupolev in 1979. Ph. D. (Eng.), deputy general director of JSC GosMKB "Raduga". Author of more than 50 publications in the field of mechanics of airspace structures.
