CC BY
75
Вычислительные технологии
Том 18, № 4, 2013
Численное моделирование течений вязкого
>к
теплопроводного газа в канале*
В. В. Шлйдуров1'2, Г. И. ЩЕПАНОВСКАЯ1, М. В. Якувович1 1 Институт вычислительного моделирования СО РАН, Красноярск, Россия 2 Университет Бейхан, Пекин, Китай e-mail: shaidurov04@gmail.com, gi@icm.krasn.ru
Предложен алгоритм численного решения уравнений Навье — Стокса для двумерного движения вязкого теплопроводного газа. Дискретизация уравнений проводится комбинацией метода траекторий для субстанциональной производной и метода конечных элементов с кусочно-билинейными базисными функциями для остальных слагаемых. Представлены результаты численного исследования структуры сверхзвукового течения в плоском канале в зоне его расширения уступом для широкого диапазона чисел Маха и Рейнольдса. Исследованы поля скоростей и давления, изучена вихревая структура циркуляционного течения в области за уступом.
Ключевые слова: уравнения Навье — Стокса, вязкий теплопроводный газ, численное моделирование, метод траекторий, метод конечных элементов.
Введение
Течение жидкости в каналах со скачкообразным расширением сечения встречается во многих технических устройствах и сооружениях. Резкое расширение сечения способно вызвать отрыв потока и существенно изменить его кинематическую структуру. Течение в плоском канале со скачкообразным расширением относится к наиболее простому классу отрывных течений c фиксированной точкой отрыва. Первые расчёты стационарных двумерных ламинарных отрывных течений несжимаемой жидкости в каналах аналитически были получены еще в 1910 г. Блазиусом в виде рядов. В дальнейшем эти расчёты использовались многими исследователями для изучения механизмов отрывных течений и тестирования разностных схем решения уравнений Навье — Стокса. В силу большой практической значимости такие течения изучались теоретически и экспериментально как для ламинарных, так и для турбулентных режимов движения несжимаемой и сжимаемой жидкости. Во многих работах данного направления рассматриваются течения в каналах с "обратным уступом" [1, 2] или с "внезапным расширением" [3-6]. Экспериментальные данные для этих случаев в плоском канале получены в [1, 3-5], где отмечается образование циркуляционной зоны за уступом. Рядом исследователей для расчётов течений с внезапным расширением были применены уравнения движения в приближении пограничного слоя.
В настоящее время ясно, что при постановке задач расчёта отрывных течений с вихревыми образованиями необходимо использовать не приближённые уравнения пограничного слоя, а полные уравнения Навье — Стокса. Однако применение более сложных
* Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00224а) и программы фундаментальных исследований Президиума РАН (проект № 18.2).
математических моделей приводит к росту вычислительных затрат. Таким примером являются задачи аэродинамики, которые на начальных этапах исследовались с помощью системы уравнений Эйлера и других приближений, а затем с развитием вычислительной техники — с использованием полной системы уравнений Навье — Стокса. Численное решение уравнений Навье — Стокса и сегодня представляет большие трудности, что обусловлено нелинейностью исходных уравнений, наличием областей больших градиентов и другими особенностями, возникающими при определённых параметрах и режимах газодинамических течений. Как следствие, возникает необходимость создания специальных численных методов решения этих уравнений. Несмотря на то что к настоящему времени разработано много численных алгоритмов и специальных комплексов программ (см. публикации [7-13] и обширную цитируемую в них литературу), проблема создания и применения эффективных численных методов и алгоритмов остаётся актуальной.
Следует отметить, что система двумерных уравнений Навье — Стокса для вязкого теплопроводного газа включает четыре дифференциальных уравнения в частных производных, вытекающих из законов сохранения массы, количества движения и внутренней энергии газа. Предложенная в настоящей работе замена искомых функций в уравнениях неразрывности и внутренней энергии переводит закон сохранения массы и полной энергии из терминов пространства Ь\ в термины гильбертова пространства Ь2. Впоследствии это значительно упрощает обоснование устойчивости и сходимости [14].
В работе для аппроксимации полной (субстанциональной, или лагранжевой) производной по времени в каждом уравнении системы используется метод траекторий, который заключается в аппроксимации этой производной с помощью разностной производной назад по времени вдоль траектории движения частицы. Под названием метода характеристик, или полулагранжевого метода, он впервые был применён в [15] для уравнения переноса массы. Далее под названием модифицированного метода характеристик он неоднократно использовался для решения уравнений параболического типа (см. работу [16] и цитированную в ней литературу). Поскольку в газовой динамике под характеристиками имеются в виду совсем другие объекты, то мы применяем более подходящее название — метод траекторий. Дискретизация по пространству остальных слагаемых уравнений Навье — Стокса на каждом временном слое проводится методом конечных элементов с кусочно-билинейными базисными функциями и применением простых квадратурных формул. Для решения систем алгебраических уравнений используется метод Якоби с улучшенным начальным приближением внутри внешних итераций по нелинейности.
Как следует из тестовых расчётов [10], модификация уравнений Навье — Стокса обеспечивает повышение точности приближённого решения по сравнению с погрешностью для немодифицированных уравнений. Вместе с тем применение комбинации методов траекторий и конечных элементов позволяет построить алгоритм, довольно эффективный с вычислительной точки зрения.
1. Постановка задачи и исходные уравнения
Рассмотрим двумерное ламинарное течение газа в плоском канале с расширением в виде уступа на нижней стенке канала при сверхзвуковой скорости потока на входе. Конфигурация расчётной области представлена на рис. 1. Начало введённой системы координат находится в левом нижнем углу в точке А. Ширина канала в левом входном сечении
У В
А
А
т С, С
и
К :
г "А
1 -►
Е
Б
Рис. 1. Канал с уступом
имеет размер Н\, в правом выходном — размер кс. Высота уступа ЕЕ — Ь = кс — Н1, длина уступа составляет с. Левая и правая границы расчётной области А1В и СП считаются достаточно удалёнными от сечения С1ЕЕ, поэтому на них можно принять условия, отвечающие невозмущённому и установившемуся течению соответственно.
Для описания движения газа используем нестационарные уравнения Навье — Стокса без упрощающих предположений. При введении безразмерных величин за масштаб длины принимается ширина канала СП, за масштаб плотности — плотность в набегающем потоке за масштаб скорости — скорость потока на входе в канал иза масштаб времени — величина Нс/и^, за масштабы давления, температуры и внутренней энергии принимаются величины из условия совершенного газа.
Выпишем дифференциальные уравнения двумерного вязкого теплопроводного газа в виде безразмерных уравнений неразрывности, количества движения и внутренней энергии
¿р
..г ди ду (Ь Р дх Р ду '
¿и
дР дтт
Р(Ь =
дх
(у дР + Р (И ду
дх
дтху
дх
+
дтх
ху
ду
+
дт
уу
ду
(е
„„ + р (ди + ду Р М \ дх ду
— д^ — дду +ф.
дх ду
Здесь ((•)/¿Ь — субстанциональная, или полная, производная, т.е.
¿Р дР дР дР
(Ь дЬ дх ду'
(2)
(3)
(4)
(5)
Р — плотность; и и у — проекции вектора скорости на оси х и у; Р = (7 — 1)Ре — давление; V = (7(7 — 1)М2е)ш — динамический коэффициент вязкости; е — внутренняя энергия. Компоненты тензора напряжений тхх, туу, тху, проекции теплового потока дх, ду и диссипативная функция Ф выражаются следующим образом:
2 / ди ду Тхх 3Ие Ц \ дх ду
уу
2 /\ру ди 3Ие Ц \ ду дх
ху
ту
ух
V (ди ду Ие \ду дх
1
де
Рг Ие дх
Чу
1
де
Рг Ие ду
х
Ч
х
ж ^
Ф = — Ие
2/ ди\2 2/ 2 (д! ди
3 \дх) 3 \ду ~ ~
2 /ди дv V
\дх ' ду) 3 \ дх ду у
(6)
где И,е — число Рейнольдса, Рг — число Прандтля, 7 = 1.4.
Для завершения постановки задачи зададим начальные и краевые условия. Пусть газ начинает движение слева направо из состояния покоя внутри области, так что р(0,х,у) = 1, и(0,х, у) = 0, г>(0,х,у) = 0. Внутренняя энергия для совершенного га-
за равна е(0, х, у) = (7 (7 — 1) М2) \ На входной границе А^Б на временном интервале
£ € (0,^п) задаются следующие параметры потока: р\
А В
1,
\АЛВ
(7 (7 — 1) М2
2Л-1
\АЛВ
0. Пусть Ъ\ = Ъ/Нс — безразмерная высота уступа, тогда профиль скорости
и(£, 0, у) на границе А^Б задается в виде
( (Ъ1 + 2а — у)(у — Ъ^/а2, у € (ЪЬЪ1 + а], и(£, 0, у) = I 1, у € (Ъ1 + а, 1 — а), (7)
[ (1 — 2а — у)(у — 1)/а2, у € [1 — а, 1),
где а — свободный параметр, который в последующих расчётах принимался равным 0.1. Выбранный профиль предназначен для обеспечения непрерывности функции и(£, х,у) в точках А1 и Б .В противном случае не только отсутствует сходимость, но и проявляются паразитические осцилляции за счёт разностного дифференцирования по пространству в окрестностях этих точек. Что касается скачка между нулевыми начальными условиями и значениями в (7) при £ > 0, то используемая монотонная аппроксимация производной по времени приводит к быстрому сглаживанию разрыва со временем.
На неподвижных твёрдых стенках выполняется условие прилипания и\г3 = 0, =0, а также условие тепловой изоляции, т. е. равенство нулю производной от внутренней энергии по нормали к твёрдой стенке де/дп\Гв, где Г = Г1 и Г3 и Г4 и Г5 — твёрдая граница. На выходе из канала в сечении СД для функций и, V, е принимаются нулевые условия Неймана, для р нет необходимости ставить здесь дополнительные условия.
е
V
2. Редукция исходных уравнений
Преобразуем уравнения (1) и (4) к новому виду. Для этого, учитывая неотрицательность плотности и внутренней энергии, введём функции
р = а2, (8)
е = е2. (9)
Произведём замену (8) в уравнении неразрывности (1) и после сокращения на 2а получим
(а 1 /ди дv\ ^
-г + - а 7— + ^ =0. (10)
(И 2 Vдх ду/ у '
Аналогично поступим в отношении уравнения внутренней энергии (4): произведём замену (9) в (4) и после сокращения на 2е получим
(е , 1 ддж 1 дqy Р (ди 1 д^ 1
Р(£ + 2е дх + 2е ду 2Д дх + ду у) + 2еФ. (11)
Используем (9) также в выражениях для теплового потока чх, чу из (6) и возьмём производные по х и у; в итоге получим
27 де
Vе'
Рг Ие дх
Чу
27 де
це-
Рг Ие ду
:12)
дЧх = 27 дх Рг Ие
де 2 д де
дх ) дх\дх)
13)
дЧу
ду
21
де
д ( де"
РгИе \Ц[ду) +еду\Цду/
:14)
С учетом (13), (14) и выражения для диссипативной функции Ф из (6) уравнение (11) принимает вид
Р
(е 7 / V (де \2 д ( де^\ 1 / °"х2
(Ь РгИе I е \дх) дх \ дх) ) РгИе \ е \ду) ' ду ду)
V (де\2 д ( де\
: к VI,-.
Р / ди ду\ 1 V 2е \ дх ду ) 2Ие е
К2+2+ + 2+2_(ди ду^2
3 дх 3 ду
дх ду 3 дх ду
:15)
Замечание. В рассматриваемой задаче внутренняя энергия положительна и больше единицы по отношению к её невозмущённой величине. Поэтому множитель 1/е не может вызвать сингулярность е "вблизи нуля" и "гасит" возможный рост давления как е2. Для совершенного газа, как следует из формулы Сазерленда, динамический коэффициент вязкости является степенной функцией от внутренней энергии, в силу чего аналогичные рассуждения справедливы для v/е.
Итак, далее будем решать систему уравнений, преобразованную к следующему виду:
(а 1 (ди ду\
(Ь +2 а[дх + ду} = 0'
16)
(и
р(Ь =
(у
р(Ь =
дР дтхх дх дх
+
дт
ху
ду
дР дт,
ду
ху
дх
+
дт
уу
ду
:17)
18)
(е 7 /V (де\2 д ( де^\ 7 /V (де\2 д ( де^ Р(Ь РгИе I е \дх) дх \ дх) ) РгИе I е \ду) ду \ду/
Р (ди ду\ 1 V 2е \ дх ду ) 2Ие е
2(2 1(дЛ2 (ду 2 2(ди дЛ2
3 \дх) 3 \ду) \дх ду) 3 \дх ду)
19)
Ч
X
2
Эту систему замыкают алгебраические соотношения для давления и динамического коэффициента вязкости совершенного газа Р = Р(а,е), V = V(а'е).
3. Метод траекторий
В качестве области определения задачи рассмотрим многоугольник П, ограниченный замкнутой ломаной ВСДЕ^А^В с границей Г, состоящей из шести сегментов:
Г1 ={ х у) х € (0.0, 10.0] у = 1.0};
Г2 ={ х у) х = 10.0, у € (0.0, 1.0)};
Гз ={ х у) х € [1.0,10.0], у = 0.0};
Г4 ={ х у) х = 1.0, у € (0.0, 0.25)};
Г5 ={ х у) х € (0.0,1.0], у = 0.25};
Гб ={ х у) х = 0.0, у € [0.25, 1.0]}.
В целях упрощения изложения примем равномерную квадратную сетку по пространству с координатами х = ¿к, у/ = ^'к, г = 0,1,...,п, ] = 0,1,...,^, и шагом к = 1/^, целиком укладывающимся по горизонтали и вертикали многоугольника П. Введённая сетка разбивает расчётную область П на квадратные ячейки ш^/ — (х, Хг+1) х (yj ,у^+1). Обозначим множество узлов этой сетки в прямоугольнике ВСДА через Бн = {si,j = (xi,yj) : г = 0,1, ...,п, = 0,1, ...,п1}, и введём сеточную область Пн = Бн П П. Обозначим через Пн = Бн П (П и Г2) множество "расчётных узлов", а через Г^ = Пн П (Г\Г2) — множество граничных узлов "известных значений" для компонент скорости. Обозначим также два участка сеточной границы как = Пн П Г2 и Г™ = Пн П Г6.
Для аппроксимации субстанциональной производной по времени в каждом уравнении системы (16)—(19) используем метод траекторий, который заключается в аппроксимации данной производной с помощью разностной производной назад по времени вдоль траектории, обусловленной уравнением (1) [13]. Для этого введём равномерную сетку по времени с шагом т = ¿дп/т:
ш7
{£к : ¿к = кт, к = 0,...,т}.
Для произвольной функции х, у) будем использовать обозначения (х, у) =
Фк,х,У) и ^ = ^(¿к,xi,Уj).
Итак, субстанциональную производную в уравнении (16) заменим разностной производной с первым порядком аппроксимации [15]:
1 - ак (X ,ук)
(20)
¿к + 1
т
где Хк = х(£к), У)к = у (¿к) — координаты траектории в момент времени £ = £к, которая при £ = ¿к+1 проходит через узел (х^у/). В принципе для определения (Хк , У)к) необходимо решить обратно по времени следующую задачу об этой траектории на отрезке £ € [¿к, ¿к+1]:
= и((>х«'уМ>' ( х((к+1 ) = Xi,
I = «(¿.хм.уй), Ьу(£к+1) =»•
Вместо этого реализуем один шаг по времени явной схемы Эйлера (тоже первого порядка аппроксимации). В итоге получим приближённые значения
Хк « Хк = ^ - тмк+1, у/ « Ук = Уj - т1>.
i,j
Ясно, что в общем случае координаты Хкк, ук не попадут в узел сетки. Поэтому значение ак(Хк,Ук) определим путём линейной интерполяции:
ак (Хк ,ук ) = а/ +
ак(х,у/) - а
к/ / Vк
х_х'
(Хк - х) +
ак(х',У) - ак.
у - у/
(Ук - у/ )
а*. - тмк-" (х,у/) - - ™к-
ак(х',У) - ак
(21)
у - у/
Координаты х и у выбираются из соображений монотонности разностной аппроксима-
ции:
х
хк-1, если и/ > 0,
у
у/-1, если V/ > 0,
х'+1 иначе, I у/+1 иначе.
В итоге монотонность (неположительность внедиагональных элементов) выполняется при условии
т < к/1 + |) для всех узлов Пн = Бн П П.
(22)
Итак, субстанциональные производные в уравнениях (16)-(19) аппроксимируются следующим образом:
(а
Р
(и
_к+1 к/ Vк Vк\
(V.
Р^
(е.
Р/
¿к + 1
¿к + 1
¿к + 1
т
и'+ 1 - ик (Хк ,Ук)
к+1 _к/
к/
Р
к+1 к/
1 - (Хк ,ук)
ек,+1 - ек (Хк ,Ук)
Р
к+1 _к/ к/
(23)
(24)
(25)
(26)
т
Значения ик (Хк ,У/к), (Хк ,У/к) и ек (Хк ,У/к) вычисляются путём линейной интерполяции аналогично формуле (21) для ак(Хк, ук).
к
4. Метод конечных элементов
После аппроксимации субстанциональной производной на каждом временном шаге £ = ¿к+1, к = 0,..., т - 1, в П и Г2 получаются уравнения
а 1 /ди .
т + 2аи + Щ) = Л' (27)
ри = - дР + д^ + д^ху + / (28)
т дх дх ду '
р^ = + ^ + ^ + /з, (29) т ду дх ду '
Ре 7 (V (де\2 д / де\ \ 7 т РгИ,е \ е\дх) дх\дх)) РгИ,е
V (де\ д / де'
е ду \ду,
¡4—
Р (ди ду\ 1 V 2е \дх ду) 2И,е е
2( ду\2 2( ду\2 /ду ду\2 2( ду ду 3 \дх) 3 \ду) \дх ду) 3 \дх ду
(30)
с правыми частями ¡1, ¡2, ¡3куда вынесены слагаемые, известные с предыдущего временного слоя.
В предыдущем разделе в явной форме были построены аппроксимации только для узловых точек, а не для всей области П и Г2. На самом деле, в методе конечных элементов после использования квадратурных формул аппроксимации в других точках не потребуется.
Для каждого узла Е Пь введём базисную функцию
(х,у)
(1 — \х — х|/^) (1 — — у\1К) ' (х,у) Е ([хг-1'хг+1] X [у^-1, у^+1]) П П'
0
иначе,
(31)
которая равна единице в в^ и нулю во всех остальных узлах Пь. Будем искать приближённое решение в виде
аН(Ь'х'у)
уН(Ь' х' у)
X] (х'у)'
АЬ'х'у)= ^ (х'у)'
£
Щ,3 (Ь)<£г,] (х'у)' еН(Ь'х'у)= X ^ (х'У)■
Для аь' иь и уь известны значения на Г^, а для еь известны лишь значения на Г™. Это следует из того, что (естественные) краевые условия Неймана в отличие от (главных) условий Дирихле не ликвидируют степени свободы в соответствующих узлах границы [17].
После стандартного использования метода конечных элементов (Бубнова — Галер-кина) с тестовыми функциями вида (31) применим квадратурную формулу трапеций для вычисления интегралов на отрезках, а её декартово произведение — для вычисления интегралов на ячейках = (х^х^) X (уз^-+1). В результате во внутренних узлах расчётной области Бь П П получим следующий сеточный аналог уравнения неразрывности (далее в этом и в следующем разделе у всех функций опущен верхний индекс к + 1' характеризующий зависимость от времени):
а
1
1
-у + (иг+1,з — и%-1,]) + -^а^з^¿¿+1 — у«-1) = ¡1 ^ е Бь П П. (32)
В приграничных или граничных узлах это и последующие сеточные уравнения для иь' уь' еь упрощаются за счёт краевых условий или меньшего носителя тестовых функций, например, на границе ГЬ"1. В итоге получим большое разнообразие уравнений, из которого ограничимся уравнениями во внутренних узлах, дающими достаточное представление о виде получаемых сеточных уравнений.
Аналогично для и приведём сеточный аналог уравнения количества движения (28) тоже лишь во внутренних узлах Бь П П:
2
+ зкк ((и"/ - Ик— 1,/)(^к—1/ + ^к/) - (и"+1/ - ик/Х^к/ + ^"+1/)) +
+ бк2Ке ((^'+1,/+1 - ^+1/— 1)^к+1,/ - (^к—1,/+1 - ^к—1,/—1)^к—1,/) + + 2к1Ке ((ик/ - и"/ —1)(Мк,/— 1 + ) - (и"/+1 - и",/Х^к/ + ^"/+1)) + + 4Л2Ке ((^'+1,/—1 - ^к—1,/—1)^к,/—1 - (^к+1,/+1 - ^к—1,/+1)^к,/+1) =
= /2 - 2к(Р'+1,/ - Р'—1/) VSi,j € Бн П П. (33)
В узлах сеточной границы в силу условия Неймана данные уравнения несколько
тлО -
упрощаются, а сеточные краевые условия на Г н вытекают из краевых условий исходной дифференциальной задачи.
Сеточный аналог уравнения (29) для компоненты скорости V во внутренних узлах Бн П П получается так же:
Рк/_ к/ + 4^1^ ((и'—1,/+1 - И'—1,/—1)^к—1,/ - (и"+1/+1 - и'+1,/—1)^к+1,/) +
+ 2к1Ке ((vi,j - Vi—1,/)(^к—1,/ + ) - - Vi,jХ^",/ + )) + 2
+ 3к2Ке ((vi,j - Vi,j—1)(Мк,/— 1 + ) - (vi,j+1 - Ч/Х^к/ + ^к,/+1)) + + ((и'+1,/+1 - ик—1,/+1)^к,/+1 - (ик+1/—1 - ик—1,/—1)^к,/—1) =
= /з - 2к (Р',/+1 - Р',/—1) V ^ € Бн П П. (34)
В узлах сеточной границы Г^ в силу условия Неймана эти уравнения выглядят проще, а значения V на Г^ являются нулевыми на основании краевых условий исходной дифференциальной задачи.
Сеточный аналог уравнения энергии (30) во внутренних узлах Бн П П получается в виде
/ / Г"/ ((е"+1/ - Гк/)2 + (ек/ - Гк— 1,/)2) -
•■к,/
т 2к2РгИ,е ек/
7 /7,- _ \2 , _ \2
2к2РгЯе ек/
Гк/+1 Гк/) + (ек/ Г",/— 1) ) +
Т
+ 2к2Р-гКё ((е",/ - Г"— 1,/)(^к—1,/ + ^к/) - (е"+1,/ - Гк/)(^к,/ + ^"+1,/)) + т
+ 2к2РгКе ((е",/ - Г",/ —1)(^",/—1 + ) - (е",/+1 - Г",/)(^к,/ + ^к,/+1)) = 1 Р" •
= /4 - (И"+1,/ - И"—1,/ + Vi,j+1 - Vi,j—1) +
+6^]^ ^ ((и"+1,/- и",/)2 +(и",/- и"—1,/)2) +
1 / / , \ 2 / \2) + 6к^ Г" ^(Vi,j+1 - ^) + - vi,j—1) ) +
1 VI 3 Г/ \2 / \ 2
+ 8№е — [(уг+1,3 — уг,3 + иг,3+1 — иг,3) + (уг,3 — уг-1,3 + иг,3+1 — иг,3) +
-г,3
+ (уг+1,з — у%,] + и^ — иг,з-1)2 + (у^з — у^-и + — иг,з-1)2} +
1 VI 3 Г/ \2 / \ 2
+ 12№е^ |_(иг+1,3 — игз — уг,з+1 + уг,3) + (иг,з — Щ-ы — уг,3+1 + уг,3) +
-г,з
+ (иг+1,з — игз — угз + угз-1)2 + (игз — иг-13 — угз + уг,3-1)2] Е Бь П П. (35)
Понятно, что в узлах ПьП(Г5 и Г2) в силу нулевого условия Неймана сеточное уравнение упрощается.
Эти уравнения дополняются краевыми условиями на Г™ из условия невозмущённого потока на входе.
5. Системы алгебраических уравнений
Для решения систем алгебраических уравнений используется многосеточный метод с внешними итерациями по нелинейности. После аппроксимации субстанциональной производной и применения метода конечных элементов с квадратурными формулами на каждом временном слое Ь = Ьк+]_' к = 0' 1'...'Ш — 1' получим системы нелинейных алгебраических уравнений.
Для сеточного аналога уравнения переноса массы система алгебраических уравнений относительно неизвестных ак+1 имеет диагональный вид
<аг,3 = Fгp3 Увз Е Пь' (36)
где ар3 — коэффициенты перед неизвестными, возможно, зависящие от других сеточных функций, а Гр3 — значения, известные с предыдущего временного слоя или из краевых условий. После исключения известных краевых значений уравнение (36) в матричном виде принимает вид
Анаь = рЛ (37)
ьь
где Аь — диагональная матрица, &ь и — векторы неизвестных и правой части размерности множества Пь.
Сеточные аналоги уравнений количества движения в терминах неизвестных значений ик+1 и ук+1 в узлах сетки Пь могут быть представлены в следующем виде:
аи иг-1,з + аи иг,з-1 + аи иг,з + аи иг,3+1 + аи иг+1,3 +
+вГ1,3-1уг-и-1 + вГ1,3+1уг-1,з+1 + [Зги+1,3-1уг+1,з-1 + = ^ ' (38)
а1'1,3 уг-1,3 + а\;3-1 угз-1 + а\:;3 у^ + а];3+1угз+1 + аг^+1,3 уг+13+ +ДГ13-1иг-1,з-1 + в~1,3+1иг-1,3+1 + ^''Ч+и- + в^+1,3+1 Щ+1,3+1 = Г^3 , (39)
\/8г.3 Е Пь.
Здесь а и в с разными индексами — коэффициенты при неизвестных и и у, равные нулю вне сеточной области Пь и зависящие нелинейным образом от других неизвестных систем (37), (41), (42) (см. ниже).
Сеточная система нелинейных алгебраических уравнений, соответствующая уравнению внутренней энергии, может быть представлена в виде
а-1 j+aj-4j-i+^++ae+lj^ = Fj v^- e пдгп, (40)
где а с разными индексами — коэффициенты при неизвестных, равные нулю вне сеточной области H и нелинейным образом зависящие от других неизвестных систем (37),
(41), (42).
Исключим известные краевые значения и последовательно занумеруем узлы H в лексикографическом порядке. Соответственно пронумеруются уравнения и неизвестные. В итоге система (33)-(35) примет матричный вид:
AUuh + = FU, + Ahvh = Fh, (41)
= Fh, (42)
где AU, A^h, Ah — пятидиагональные, B^B", — четырехдиагональные матрицы; uh, vh, eh и FU, F^h, Fh — векторы неизвестных и правых частей размерности множества П. Особенностью как системы (41)-(42), так и её промежуточных линеаризаций является необходимость одновременного использования трёх матричных уравнений во всех видах итераций (итераций Якоби для линеаризованных систем и итераций по нелинейности).
Таким образом, получена вариационно-разностная схема первого порядка аппроксимации по времени и по пространству. Для решения систем линейных алгебраических уравнений на каждом временном слое применялся точечный метод Якоби [18]. Сходимость этого метода и итераций по нелинейности значительно ускоряется при использовании в качестве начального приближения квадратичной экстраполяции значений по времени с двух временных слоев вместо простого переноса значений с предыдущего слоя. Ввиду существенного диагонального преобладания среднее количество итераций, необходимое для сходимости метода Якоби на сетке 1001x101 узлов, составляло не более 10.
6. Расчёт течения газа в канале с различными числами Маха и Рейнольдса
Приведённый алгоритм реализован для сформулированной выше задачи течения газа при сверхзвуковой скорости на входе. В качестве уравнений (5) в расчётах использованы уравнения состояния совершенного газа
P = (Y - 1)pe, T = y(y - 1)M2e,
зависимость динамического коэффициента вязкости газа представляется формулой Са-зерленда
ß = Tш.
Различные модификации этих уравнений и условия их применения приведены в [13].
Расчёты выполнялись на сетке, содержащей 1001x101 узлов, шаг по пространству h = 0.01, шаг по времени т = 0.001. Газодинамическая постоянная y, число Рейнольдса Re, число Прандтля Pr, число Маха Ми ш имели следующие значения: y = 1.4, Re = 2 ■ 103 и 104, Рг = 0.72, M = 2 и 4, ш = 0.8.
Рис. 2. Продольная составляющая скорости в моменты времени £ = 4 (а) и £ = 50 (б)
0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3
Рис. 3. Изолинии плотности для моментов времени £ = 4 (а), £ = 50 (б); I — И = 2, Ив = 2 • 103, II - И = 4, Ив = 104
На рис. 2 показана картина течения в канале для М = 2, Ие = 2 • 103. На рис. 2, а представлена компонента скорости и в момент времени Ь = 4, на рис. 2, б — в момент времени Ь = 50. За уступом происходит формирование вихря с отрицательными значениями скорости. С течением времени вихревая зона увеличивается в направлении потока и за уступом формируется течение со скоростью, близкой к нулевым значениям. Следует отметить, что за характерный размер Ь принята ширина канала. В данном случае, чтобы точка "примыкания" основного потока к донному течению оставалась в пределах расчётной области [19], длина канала при расчётах равнялась 10Ь.
На рис. 3 приведены соответствующие изолинии плотности для тех же моментов времени. Видно, что для Ь = 4 характерны формирование отрывной зоны в районе уступа и присутствие головного скачка. Со временем распределение плотности демонстрирует установление течения.
Таким образом, замена искомых функций в уравнениях сохранения массы и внутренней энергии привела к меньшей абсолютной погрешности в нормах Ь2 и Ьж, что и раньше было отмечено в одномерном случае [10]. Применение комбинации метода траекторий и метода конечных элементов не требует согласования триангуляций на соседних временных слоях. Это значительно облегчает динамическое разрежение или сгущение
триангуляций по времени для оптимизации вычислительной работы или улучшения аппроксимации в пограничных слоях и ударных волнах. Для решения систем алгебраических уравнений ввиду значительного диагонального преобладания использовался метод Якоби в комбинации с внешними итерациями по нелинейности. Совместное использование методов траекторий и конечных элементов позволило построить вычислительно устойчивый и экономичный алгоритм.
Список литературы
[1] Armaly В.Р., Durst F., Pereira J.C.F., Schonung В. Experimental and theoretical investigation of backward-facing step flow //J. Fluid Mech. 1983. Vol. 127. P. 473-496.
[2] Елизарова Т.Г., Соколова М.Е., Шеретов Ю.В. Квазигазодинамические уравнения и численное моделирование течений вязкого газа // Журнал вычисл. математики и матем. физики. 2005. Т. 45, № 3. С. 545-556.
[3] Cherdron W., Durst F., Whitelow J.H. Asymmetric flow and instabilities in symmetric ducts with sudden expansion //J. Fluid Mech. 1978. Vol. 84. P. 13-31.
[4] Durst F., Melling A., Whitelow J.H. Low reynolds number flow-over a plane symmetric sudden expansion // Ibid. 1974. Vol. 64. P. 111-118.
[5] Fearn R.M., Mullin T., Cliffe K.A. Nonlinear flow phenomena in a symmetric sudden expansion // Ibid. 1990. Vol. 211. P. 595-608.
[6] Hawa T., Rusak Z. The dynamics of a laminar flow in a symmetric channel with a sudden expansion // Ibid. 2001. Vol. 436. P. 283-320.
[7] Bosnyakov S., Kursakov I., Lysenkov A. et al. Computational tools for supporting the testing of civil aircraft configurations in wind tunnels // Progress in Aerospace Sci. 2008. Vol. 44. P. 67-120.
[8] Ковеня В.М., Слюняев А.Ю. Моделирование сверхзвуковых течений газа в канале // Вычисл. технологии. 2007. Т. 12. Спец. выпуск. 4. С. 32-41.
[9] Oberkampf W.L., Trucano T.G. Verification and validation in computational fluid dynamics // Progress in Aerospace Sci. 2002. Vol. 38. P. 209-272.
[10] Шайдуров В.В., Щепановская Г.И., Якубович М.В. Применение метода траекторий и метода конечных элементов в моделировании движения вязкого теплопроводного газа // Вычисл. методы и программирование. 2011. Т. 12. С. 275-281.
[11] Vos J.B., Rizzi A., Darracq D., Hirschel E.H. Navier — Stokes solvers in European aircraft design // Progress in Aerospace Sci. 2002. Vol. 38. P. 601-697.
[12] ADIGMA — A European Initiative on the Development of Adaptive Higher-Order Variational Methods for Aerospace Appl. Vol. 113. Notes on Numerical Fluid Mechanics and Multidisciplinary Design. Springer, 2010. P. 339-353.
[13] The Handbook of Fluid Dynamics / Ed. R.W. Johnson. CRC Press LLC & Springer, 1998.
[14] Ушакова О.А., Шайдуров В.В., Щепановская Г.И. Метод конечных элементов для уравнений Навье — Стокса в сферической системе координат // Вестник КрасГУ. 2006. № 4. 151-156.
[15] Pironneau O. On the transport-diffusion algorithm and its applications to the Navier — Stokes equations // Numer. Math. 1982. Vol. 38. P. 309-332.
[16] Chen H., Lin Q., Shaydurov V.V., Zhou J. Error estimates for trangular and tetrahedral finite elements in combination with a trajectory approximation of the first derivatives for advection-diffusion equations // Numer. Analysis and Appl. 2011. Vol. 4, No. 4. P. 345-362.
[17] Rannacher R. Methods for Numerical Flow Simulation. Institute of Appl. Math., Univ. of Heidelberg, Germany, 2007. 58 p.
[18] ВоЕводин В.В., Кузнецов Ю.Л. Матрицы и вычисления. M.: Наука, 1984.
[19] Блшкин В.А., Егоров И.В., Иванов Д.В. Торможение сверхзвукового потока в плоских и осесимметричных каналах // Изв. РАН. Механика жидкости и газа. 1998. № 2. С. 143-152.
Поступила в 'редакцию 6 марта 2013 г.
