CC BY
106
Том XLIII
УЧЕНЫЕ ЗАПИСКИ ЦАГИ 2012
№ 6
УДК 519.688 533.69
ПОВЫШЕНИЕ УСТОЙЧИВОСТИ ЯВНОЙ СХЕМЫ ГОДУНОВА — КОЛГАНА — РОДИОНОВА ЛОКАЛЬНЫМ ВВЕДЕНИЕМ НЕЯВНОГО СГЛАЖИВАТЕЛЯ
Е. В. КАЖАН
Рассмотрена возможность ускорения стационарного численного решения задач обтекания ЛА потоком вязкого газа по явной схеме Годунова — Колгана — Родионова при помощи локального введения неявного сглаживателя. Сглаживатель вводится только в наиболее жестких, с точки зрения устойчивости явной схемы, областях, например, во внутренних областях пограничного слоя, где размер ячеек наименьший.
Ключевые слова: вычислительная аэродинамика, неявная схема, схема Годунова — Кол-гана — Родионова.
ВВЕДЕНИЕ
В ЦАГИ была разработана и получила широкое применение технология вычислительного эксперимента под названием «Электронная аэродинамическая труба» (ЭАДТ) [1, 2]. Разработанный для этого пакет прикладных программ (ППП EWT — ЦАГИ) [3] включает весь набор инструментов для реализации вычислительного эксперимента. Вычислительный эксперимент позволяет исследовать турбулентные течения вязкого газа при обтекании реальных конфигураций, в том числе моделировать реальный аэродинамический эксперимент с учетом особенностей экспериментальных установок. В данной статье предложено, как в рамках данной технологии ускорить получение стационарных решений. ППП EWT — ЦАГИ включает в себя несколько солве-ров, в том числе V3Solver [4] и ZEUS [5]. Для получения стационарных решений в этих солверах используется метод установления по времени и стационарное течение получается как предел нестационарного. Все численные схемы, которые могут быть использованы в методе установления, можно разделить на два класса — явные и неявные. Они имеют свои достоинства и недостатки.
Неявные схемы могут быть абсолютно устойчивы, но требуют больше вычислительных затрат на совершение итерации, чем явные схемы. И явные, и неявные численные схемы широко используются в вычислительной аэродинамике и имеют множество реализаций [6]. В ЭАДТ до настоящего времени использовались только явные схемы. Для ускорения получения стационарных решений применялся метод локального шага по времени. Это означает, что в разных ячейках сетки расчет проводится с разными значениями шага по времени, определяемыми локальными условиями устойчивости. Существуют два подхода, которые позволяют получать стационарное решение быстрее, чем при использовании явной схемы с локальным шагом по времени. Один подход — это многосеточные методы [7, 8], другой — использование неявной схемы. В данной работе был выбран второй подход.
Современные неявные методы решения стационарных задач вычислительной аэродинамики, как правило, основаны на линеаризации
КАЖАН Егор Вячеславович
младший научный сотрудник ЦАГИ
зависимости от времени уравнений движения газа [6]. Значения параметров на неизвестном временном слое (n +1) представляются в виде:
— С)— ( )+| ( )Au
——
где Дм = un+1 — un, [un j — некоторая аппроксимация матрицы Якоби, вычисляемая по пара-
ди v '
^ n Т> ^
метрам на известном временном слое и . В результате, аппроксимация уравнений движения сводится к системе линейных алгебраических уравнений (СЛАУ) относительно приращения параметров Дм (неявная схема в «дельта-форме» [9]). При решении трехмерных уравнений Навье — Стокса и Рейнольдса матрица такой СЛАУ имеет блочно-ленточную форму с несколькими ненулевыми блочными диагоналями, разделенными многочисленными нулевыми диагоналями. Существующие методы решения таких СЛАУ можно подразделить на два больших класса.
В первом классе методов СЛАУ для Дм не упрощается. Несколько «шагов по времени» с точным решением СЛАУ при бесконечно больших шагах по времени эквивалентны решению нелинейной стационарной системы уравнений методом Ньютона, что обеспечивает квадратичную скорость сходимости при наличии хорошего начального приближения [11]. Однако, как правило, начальное поле сильно отличается от решения, поэтому итерационный процесс приходится начинать с небольших шагов по времени. При точном решении СЛАУ требуемые вычислительные ресурсы неприемлемо велики [6]. Весьма употребительны различные версии метода сопряженных градиентов [11, 12], среди которых особенно эффективным является обобщенный метод минимальных невязок GMRES [13] в безматричной реализации [14]. Однако в общем случае для таких методов решающее значение имеет применение различных предобусловливателей [15]. При решении плохо обусловленных практических задач эти методы не обеспечивают достаточную робастность [15].
Более простым в реализации является второй класс методов, основанный на упрощении СЛАУ. Здесь можно выделить методы факторизации, связанные с представлением матрицы системы уравнений в виде произведения нескольких матриц специального вида, и методы, связанные с представлением матрицы системы уравнений в виде суммы нескольких матриц. Высокоэффективная схема первой группы была предложена в работе [16]. Это так называемый метод переменных направлений (ADI), при котором матрица системы факторизуется на три множителя, связанных с численным дифференцированием в каждом из пространственных направлений. Каждый множитель приводится к системе, решаемой матричной прогонкой или даже рядом скалярных прогонок. Однако метод ADI в трехмерном случае обладает условной устойчивостью. Существует ряд других способов факторизации, обладающих абсолютной устойчивостью, например [17]. Во второй группе схем, основанных на упрощении матриц системы уравнений, используются такие эффективные методы решения СЛАУ, как блочные методы Якоби и Гаусса — Зей-деля [18] и методы релаксации, среди которых особенно популярен метод LU-SGS и его модификация LU-SSOR [19]. Как правило, число итераций таких методов на каждом шаге по времени ограничивается [6], так что решение упрощенной СЛАУ не находится, что не мешает сходимости процесса в целом к стационарному решению. В отличие от методов первой группы методы этого типа легко обобщаются на неструктурированные сетки [20]. Смена направления обхода ячеек структурированной сетки, которая ускоряет сходимость метода Гаусса — Зейделя, на неструктурированных сетках может быть успешно заменена произвольной перенумерацией ячеек [21].
При упрощении СЛАУ неявной схемы существенную роль играет сокращение шаблона схемы на неявном слое. Наиболее употребительным способом для достижения этой цели является метод отсроченной коррекции [22], при котором неявный оператор, применяемый к Дм, аппроксимируется с первым порядком точности по пространству на компактном шаблоне, и лишь в явном операторе применяется аппроксимация высокого порядка на развернутом шаблоне. Метод отсроченной коррекции основан на том, что при приближении к стационарному решению неявный оператор стремится к нулю. Это обеспечивает высокий порядок аппроксимации стацио-
нарного решения. При этом использование неявного оператора первого порядка точности повышает надежность схемы.
Предлагаемый в настоящей работе неявный метод основан на методе отсроченной коррекции и использует для решения системы линейных уравнений блочный метод Гаусса — Зейделя с перенумерацией ячеек. Этот метод продолжает традиции отечественной школы вычислительной аэродинамики, основанной С. К. Годуновым и его последователями, в частности — А. Н. Крайко и С. М. Босняковым. Ключевым принципом этой школы является использование при аппроксимации системы дифференциальных уравнений ее математических свойств, и прежде всего, учет направления распространения информации. В явной схеме Годунова [23] для этой цели впервые было использовано решение задачи Римана о распаде произвольного разрыва. Неявная схема, предложенная в классической книге [24], не относится к этому классу, равно как и базовый вариант метода ADI [16]. Однако уже в одной из последующих работ авторов ADI [17] был предложен метод ADI с учетом направления распространения информации, основанный на расщеплении вектора потоков. Родственными работе [17] являются классические отечественные работы [25, 26]. В них также используется подход ADI, матрицы Якоби конвективных потоков аппроксимируются на гранях ячеек, а в явном операторе используется нелинейная схема Годунова второго порядка аппроксимации. Линеаризация Роу [27] признана оптимальным компромиссом между точностью и эффективностью при вычислении параметров на гранях ячеек [6]. В работах [28, 29] описана неявная схема, основанная на решении задачи Римана без использования факторизации неявного оператора, с решением СЛАУ блочным методом Гаусса — Зейделя. Автор статьи опирался на опыт применения подобных схем разработчиками программы [30].
При построении неявных методов для решения уравнений Рейнольдса, замкнутых дифференциальной моделью турбулентности, возникает проблема выбора аппроксимации источников, связанных с производством и диссипацией турбулентности. Нередко для этого используется неявная аппроксимация, однако известно, что она устойчива не во всех случаях [6]. В работе [31] предложена абсолютно устойчивая схема, в которой для производства турбулентности используется явная аппроксимация, а для диссипации — неявная. В данной работе используется другой метод, основанный на анализе собственных чисел матрицы Якоби для источников. Этот метод был сформулирован автором самостоятельно как логическое развитие идей, изложенных в [32]. Впоследствии автору стало известно, что аналогичный подход был предложен в работе [33], где сообщается, что он обеспечивает более быструю сходимость, чем метод [31].
Главной оригинальной особенностью численного метода, который описан ниже, является локальный выбор явной или неявной схемы и способа осуществления шага по времени (глобальный, локальный) в зависимости от соотношения между заданным глобальным шагом по времени и локальным условием устойчивости явной схемы. Далее будет показано, что данный комбинированный подход позволяет существенно ускорить получение стационарного решения по сравнению с методами, использующими одну и ту же схему во всей расчетной области. Следует отметить, что существует обширный класс методов зональной декомпозиции [34], где в разных блоках расчетной сетки используются разные численные методы. В отличие от методов этого класса в предлагаемом подходе выбор схемы определяется в каждой ячейке сетки на основе локальных параметров численного решения, при этом удается сохранить однородность алгоритма. Локальный выбор схемы обеспечивает большую гибкость метода. Эффективность данного подхода продемонстрирована на ряде тестовых задач.
БАЗОВАЯ ЯВНАЯ СХЕМА ГОДУНОВА — КОЛГАНА — РОДИОНОВА
Решается система уравнений Рейнольдса, замкнутая моделью турбулентности q-ю [35—37]. Численная схема строится в рамках конечно-объемного подхода.
Используется маршевый метод расчета по времени. Для получения стационарного решения используется метод установления, когда задается некоторое начальное приближение, а стационарное решение получается как предел нестационарной адаптации течения к заданным стационарным граничным условиям.
Применяются многоблочные структурированные сетки с шестигранными ячейками. В ходе расчета используются параметры газа в центрах ячеек сетки. В общем виде аппроксимация законов сохранения записывается в виде
_п
а"*' = - - £ Р° + тпЖ,а. (1)
— дЕ,
Здесь Ра — потоки величины а сквозь грань дЕ,. При описании течений газа рассматриваются следующие законы сохранения: 1) закон сохранения массы (а = р); 2) закон сохранения каждой из трех компонент импульса (а = ри; ру; р>^); 3) закон сохранения энергии (а = рЕ). Модель турбулентности д-ю включает в себя два дополнительных дифференциальных уравнения для параметров турбулентности — характерной величины пульсаций скорости д и характерной частоты турбулентных пульсаций ю. В этом случае к системе законов сохранения добавляются еще два (а = рд; рю). Объединяя аппроксимации (1) всех законов сохранения в одну систему и
—► —►
вводя вектор консервативных переменных и, векторы их потоков Р сквозь грани дЕ, ячейки ,,
—►
а также вектор источниковых членов Ж, получим общую формулировку численной схемы:
ип+1 = ип -—£ Р + тпЩ . (2)
, дЕг- источники
потоки
—►
Вектор источниковых членов Ж описывает производство и диссипацию параметров турбулентности. Особенности рассматриваемой базовой явной схемы, подробно изложенной в [32]: для описания конвективных потоков используется явная монотонная схема Годунова — Кол-гана — Родионова [24, 38, 39], для описания диффузионных членов — центрально-разностная аппроксимация с модификацией для неравномерных сеток, а для источниковых членов — локально-неявная аппроксимация. Схема номинально имеет второй порядок аппроксимации по пространству и времени.
ЛИНЕАРИЗОВАННАЯ НЕЯВНАЯ СХЕМА НА ОСНОВЕ СХЕМЫ ГОДУНОВА — КОЛГАНА
—►
В уравнении (2) вектор и" известен (он вычисляется через параметры газа на известном временном слое п). Вектор и,п+1 нужно найти. Для этого необходимо указать способ вычисления
суммы потоков £ Р и источниковых членов Ж. Можно получить систему нелинейных алгеб-
дЕ,
—►
раических уравнений, трактуя соответствующие соотношения между потоками Р как нелинейные функции от известных векторов и, и искомых и,"*1. При линеаризации этих соотношений искомый вектор и,"*1 можно найти как решение СЛАУ. Если же вычислять эти потоки только по
известным величинам ип, схема будет совпадать с явной схемой Годунова — Колгана [24, 38]. Рассмотрим отдельно схемы для конвективных, диффузионных и источниковых членов. Схема для решаемой системы уравнений Рейнольдса будет линейной комбинацией этих схем.
Вначале рассмотрим неявную схему с линеаризацией для упрощенной системы уравнений Рейнольдса, замкнутой моделью турбулентности д-ю, в которой оставлены только конвективные члены:
ї +$1«
V дЕ
= 0,
и =
Обозначения:
^ , у
РЕ =
Р(и
def р I и2 + V2 +
р ^ ' рУп '
ри рупи + РЯх
Pv р^ + ряу
р^ —► , ./конв = PУnW + РЯ2
рЕ Уп (рЕ + Р)
рЧ руп4
р® V V руп® V
ч2 + 1 ^p, К — 1 def Уп = иЯх + + wsz,
(3)
где
3 = (х, 5^, 52) — вектор внешней нормали к грани ячейки дЕ,, длина которого равна площа-
ди грани.
Неявная схема имеет вид
-*п+1 ->п
Уи, - и ,
І п
+
п+1
Р = 0.
дЕ,
Чтобы максимально сохранить особенности базовой явной схемы аналогично [25, 26], перепишем ее в форме приращений относительно шага по времени. Введем оператор А:
ап. Тогда:
Аап = ап+1
у, Ап+£ Рп + £ аР = 0.
(4)
дЕ,
дЕ,
При такой форме записи первые два слагаемых (4) полностью совпадают с явной схемой (2) при условии отсутствия источников и потоков молекулярной диффузии.
Вся «неявность» заключена в последнем третьем слагаемом уравнения (4). Очевидно, что при сходимости к стационарному решению первое и последнее слагаемые (4), являющиеся разницей параметров на явном и неявном временных слоях, исчезают, и стационарное решение полностью определяется аппроксимацией потоков на явном слое. Аналогично [25, 26] применим также разные аппроксимации для потоков на явном слое и для приращений потоков. Аппроксимацию для потоков из второго слагаемого выберем такую же, как в явной схеме, уже реализованной в солверах ЭАДТ, а для аппроксимации приращений потоков выберем способ, допускающий линеаризацию — аппроксимацию Роу [27]. Следует отметить, что аппроксимация потоков из второго слагаемого имеет второй порядок, а аппроксимация приращений потоков — первый.
Потоки через грань дЕ,, согласно [27], запишем в линеаризованном виде
рп+1 РдЕ, = Ап = АдЕ, Цп+1 и дЕ, + СПпЕг
рп рдЕ, А = иЕ +СдЕ,
где А-матрица Роу, С — константа. В форме приращений:
аре = аПе АЦП,
А ~*п
У Ап + £ АпАи + £Рп = 0.
Т дЕ,- дЕ
Величины на гранях ищ, Цщ и прочие находятся из решения задачи о распаде произвольного разрыва по Роу, поэтому для них можно записать соотношения процедуры Роу через параметры слева иь и справа от разрыва ик:
-*п , -*п
uRdEi + uLdEi
і
-*п -*п
uRdEi — uLdEi
2 I $Ei 2 ^ 2 1 ldE 2
В этих формулах |l| = P(sign(Л))P-1, где РЛР-1 = A; P — матрица, столбцами которой являются собственные векторы матрицы Роу A; Л — диагональная матрица собственных чисел матрицы A. Для линейности будем считать матрицу |l||E, приближением для матрицы |l^ •
—►
Тогда приращение вектора U представится в виде
n+1 - n+1
Un+1 _ URdEi + ULdEi UdEi __
I
^ n+1 ^ n+1
n+1 URdEi — ULdEi dEi
WE _
AURdEt + AuLdEi |i in AuRdEi ~AuLdEi
2 I Щ 2
Пусть индекс NBR обозначает соседнюю ячейку, имеющую общую грань dE, с ячейкой i. Чтобы получить неявную схему первого порядка, для левых граней используются значения
UhdEi = UNBR, URdE, = ui, а для правых граней ULdE, = ui, URdE, = UNBR.
Схема, учитывающая только конвективные члены, будет выглядеть следующим образом:
Rконв i Aui + £ Rконв NBR AuNBR
dEi
_—-£ Fn,
V z_i
i dEi
(5)
где введены обозначения:
1+T1
Vi
£
dE+i
An +
£
dE — i
An —
Чонв NBR AuNBR _ y dEi Vi
' An — An ' An + An
£ AuNBR — £ AuNBR
2 2
dE+i z, К J dE—г А К J
Здесь дЕ +, — грани, внешние нормали к которым направлены в сторону увеличения индексов ячеек структурированной сетки, дЕ —, — грани, внешние нормали к которым направлены в сторону уменьшения индексов ячеек структурированной сетки. |Л| = Р | Л | Р—1, | Л | — диагональная матрица, у которой на диагонали стоят модули собственных чисел матрицы Роу Л.
Пусть совершено п шагов по времени, и мы делаем (п + 1). На п слое по стандартной («явной») процедуре Годунова — Колгана второго порядка аппроксимации по пространству [24, 38] вычисляем потоки через боковые грани и их сумму £ Рп. На п-м слое вычис-
дЕ
ляются осреднения Роу скорости на гранях (и, V, ^ и энтальпии на гранях
def
e + p'
2 2 2 u 2 + v2 + w2
+—— P ) по формуле:
л/рТ+JpR
VpTvl +4pRv л/рТ + y[pR
_ л/ріУг +'ІРя™я _4рьнь + ч/рКнк
л/рТ+ ЛК
л/рТ+ ТрК
По этим значениям вычисляются компоненты матрицы Роу А на гранях.
По матрицам А однозначно вычисляются матрицы |А| = Р | Л | Р-1. Удобно вычислять
не матрицу Роу А, а сразу собственные векторы и собственные числа. Все компоненты этих матриц вычисляются на явном слое.
Для учета вязкости в потоки (3) добавим диффузионные члены:
1 йёУ+/¿г _ 0, / _ 7К0НВ+7ДИфф,
дґ
дЕ
^ т
7дифф _(° ( (п (п {Іхпй + (пУ + 11ПК + °п ) ) Тп ) ,
где 1^ — компоненты тензора диффузионных потоков импульса; си — диффузионные потоки тепловой энергии; ТЩ и Ти® — диффузионные потоки параметров турбулентности [36].
Для Iлинеаризация даст приближение:
1ш1 _ С + АІххгх + АІхугу + ЫХ1гг, Іхх _ | 2 рд2 | - ц
4 дй 2 Эу 2 дw 3 Эх 3 Эу 3 дг
Іху _ —
дй + ду ду дх
, Іхг _-
дй + дw дг дх
М М-лам + Мтурб •
Рассмотрим производные неконсервативных переменных д//дх в криволинейной системе координат:
д/ Ж К+7 дп+7 К дх Э^ дх Эп дх Э£ дх
Проведем упрощение [40] — пренебрежем производными по криволинейным координатам, касательным к граням:
( 7 I д/ э^ /ИВК - / э% ,
Эх, у, г )дЕ Э£, Эх, у, г А^ВК + А, дх, у, г ’
где Лг — расстояние между центром ячейки и центром грани дЕ,. Итого:
(
АІххдЕ, _-МпЕ
4 Э^ Э^ Э^
---------гх +---------+--------------гу
3 Эх ду у дг
1 АйЫВК -Айі
ЭЕ,- А^ВК + Аг
АІхудЕ, _-МпЕ,
2 I Ау«вк -Ау,
------------гх +-----------
3 Эу дх
дЕ, Аывк + А ,
п , 2 Э^ Э^
АІхгЭЕ, _-МпЕ, | -Т^Т~ гх + ^~ г
АwЖК - А^'
3 д дх )дЕ, А 1Нвк +А,
Аналогично вычисляются производные по другим координатам для членов из уравнений у и г-компонент импульса для уравнения энергии.
В уравнение энергии входит слагаемое вида
(хпи + ( + Ігп™ + (п + Р°п ) ).
Линеаризация дает приближение:
ЭЕ +1хп ПЕі АидЕ,
(хпи Т+1 =(хпи Т. +1
ЛиДЕ. + Д!хп дЕ.и "е. ; (п + Р°п -(п + Р°п Ур. + Д(п + Р°п),
ЭЕ.
Члены вида 1хп щ ДмдЕ. в приращении потока энергии связаны с работой вязких и турбулентных напряжений при конвективном переносе газа. Поэтому величины ДмЭ£., Дуд£., Д^е следует вычислять так же, как и при вычислении конвективных потоков, т. е. по распаду Роу. —► —► —►
Обозначим тогда ДЕЦифф - Д^Еида + ДЕда.
Др
Еит
- (о 0 0 0 (п (ЛидЕ + Іуп Щ ЛуЭЕ. + Ігп Щ Л™ЭЕ.
о о
-
Л1х
ЛІ
уп
ЛІг
(хпи + (уп^ + ДІгп™ + ■Л(Ч п + Р°п ))
лтщ
■С
Применяем линеаризацию и выражаем значения вектора консервативных величин на гранях
„ Е + II _ Е - II _
ЛидЕ -----ЛиКЭЕ +--------Ли^дЕ- через параметры слева и справа (индексы Ь и К). Тогда:
ЛР
ЕиШЭЕ,
ЭР
ЕиЖ
Эи
ЛиЭЕ. - (%+ )^Е. ЛиЬдЕі + (- Т^)Е. ЛиКдЕі,
ЭЕ.
г±-
Э^ЕШ8 Е ± И = Э^ЕШ8 Е ± Р ( (ЛТ)Р-1
Эи
Эи
ЭР
ЕШТ8
Эи
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0
Си+ І> + 4> Іп хп I” уп г гп 0 0 0
Р Р Р Р
0 0 0 0 0 0 0
0 0 0 0 0 0 0
Для Л^£ а#, осталось расписать приращения потоков диффузии тепла:
( .. п \
Ä(n + Рсп )-. =-
к
к — 1
Млам dEj + Мтурб dEj
Рг!м Prn
турб
,Э4„ ,Э4„ ^ ät^r — ÄTj.
—sx +------sy +-s
v dx dy dz ' у
dE, ÄNBR + Äi
ÄTnq = —
Mn + Мтурб ЭЕ,
Млам dEi +
V
Prq ,ro турб у
ÄT =■
L-s~Ln
Mn + Мтурб dEi Млам dE, +
V
Prq,ro турб у
Э4 Э4 Э4
—sx + —sy + —sz dx dy y dz
Э4 Э4 Э4
—sx +----sy +-----sz
Эх x dy y 3z z
ÄqNBR — Äq j
dE, ÄNBR + Äi
ÄroNBR — Aro,
dE, ÄNBR + Äj
Теперь приращения диффузионных потоков через грань dEi можно выразить через приращения консервативных параметров в ячейках:
AFNSdE, = AdEi (^)(rNSßAuNBR - ГAui )•
Здесь
Г Э(и v w T q ю —
Эи
— матрица преобразования от консервативных пере-
менных к неконсервативным, а А(^) — матрица коэффициентов, зависящих от криволинейных координат. Заметим, что матрица преобразования Г. не квадратная — так же, как и матрица А (^). Квадратной матрицей будет их произведение.
Величины с явного слоя вида цПе- интерполируются на грань по параметрам в соответствующих ячейках. Параметры рПЕ. , ыПе , v¡¡E , WnE берутся из величин на грани с явного временного слоя. Величины на грани вычисляются из решения задачи Римана по Роу:
^ П ^ П
иЯ дЕ. - иЬ дЕ.
Un = UR dEj + UL dEj |T| n
U dEj =---------------------2---------------І1 dEj
-, где, по-прежнему, IT = P (sign(A) )P 1.
Таким образом, чтобы учесть молекулярную вязкость, достаточно в (5) добавить к матрицам Я поправку:
R = R + R
конв дифф ’
R
дифф i
= £[(( — Z—nE—, )-(+, (4) — AЖ—, (4)-
dEj
Аналогично, требуется поправить и матрицы R nbr :
£ Rдифф NBR = ( +VE — v AdE—, (4)ГNBR ) + (Z—nE + + AdE+, (4-Г
NBR
dE,
Рассмотрим аппроксимацию источников.
При использовании модели турбулентности д-ю в (3) войдут следующие источники:
dt J UdV + <j) fds =(0 0 0 0 0 Wq Wro-,
dE
Wq =pq
(Aq Cq pro ^ Cro(pro—
p + Bq +-^ ™ -~A ' D — ' roVh^ ’
pro p
, Wro pAro + Bropro +
Аппроксимация и линеаризация дают приближение:
Жп+1 = Жп +
д д 1
¿К джп
- Ар- + - Ар сс,
Эр<?
Эрос
дЩп
ЩП+1 = Ж” + ^с Ар сс.
ю
Эрос
Матрица Якоби источников имеет вид:
(0 0 0 0 0 0
0 0 0 0 0 0
Э\¥
дії
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 х„
ч
0 0 0 0 0 0
( а
СО
0
0
0
0
0
-+с
2 -
к
дЩ А- дЩ
где х-=^—=—+в- + С с
4 др- (й дрос
= Вс + 2Сс сс — собственные числа этой матрицы.
Источники порождают экспоненциально растущие моды решения вида е^, где X, — собственные числа данной матрицы. Явная аппрокисмация абсолютно устойчива при описании растущих мод, но условно устойчива при описании затухающих мод. Наоборот, неявная аппроксимация абсолютно устойчива при описании затухающих мод, но условно устойчива при описании растущих мод [41].
«Характеристическими» переменными для источников уравнения (6) будут: р, ри, ру, pw, рЕ,
( А ^
~+с
р-
со
Х Ю-Ч
-р сс, р сс. Введем аппроксимацию источников с параметром :
И-+1 = Ж-п + Х П и-А
р---
( А >
—-+с ^ -со
Л п Л п
Х ю Х-
Л ( ^ —-+с
р
+ -
со
т, п т, п
Х с Х-
=+х П и-Ар-+
( А,
со
- + С
2 -
X си т - X П и,,
п
^ іп 1П
Х с — Х -
Ар сс,
щп+1 = < + Х с и „Ар сс.
1, хп < 0
0, х\ > 0,
иш=‘
1, хс < 0 0, хп > 0
дают абсолютно устойчивую аппроксимацию,
зависящую от знака собственных чисел. Аналогичная аппроксимация была предложена в [33]. Если же взять и- = 0, и щ = 0, получим явную аппроксимацию, абсолютно устойчивую при
положительных собственных числах. Для получения абсолютно устойчивой аппроксимации нужно добавить в (5) к матрице Кконв г- следующую матрицу:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 -Хпаиа
ч ч
0 0 0 0 0 0
—-2+с
(wn )2 -
n Ä.X-Ä,- и-
Л У1 Л П
Хю - X -
-XX
И явная, и абсолютно устойчивые аппроксимации дают в вектор правой части добавку:
W = (0 0 0 0 0 X- (р-)n pn (-C”(ron)2 )(рю)
ВВЕДЕНИЕ НЕЯВНОГО СГЛАЖИВАТЕЛЯ
Для системы уравнений Навье — Стокса, замкнутой моделью д-ю турбулентности, схема, аналогичная (3), представима в виде:
R Au, +I R NBR Au
ЭЕ,
NBR =-VI Fn +TnWn,
V ЭЕ
(7)
где для каждой ячеики введены квадратные плотные матрицы, по размерности совпадающие с количеством неизвестных в ячейке:
R; -^конв і + -^дифф і + Rsource і , RNBR ^KraiiENeR + -^дифф NBR •
конвекция диффузия источники
конвекция диффузия
Следует помнить, что каждая строчка (7) на самом деле содержит векторы консервативных переменных и их произведения на матрицы, т. е. представляет собой блок из нескольких строчек матрицы системы линейных алгебраических уравнений.
Итерационные методы решения данной СЛАУ, рассматривающие итерационные значения приращений консервативных параметров в соседних ячейках как известные, позволяют сохранить структуру явного метода. Явную схему формально также можно записать в виде системы линейных алгебраических уравнений (правда, матрица этой системы будет единичной). И для решения этой системы также можно применить итерационный метод, дающий, правда, абсолютно точное решение системы за одну итерацию. Обозначим верхним индексом k номер итерации. При этом станут наглядными все отличия явной схемы от рассматриваемой неявной:
R,Auf =-I RNBRAUNBr"-6 -—£ Fn +TnWn — неявная схема,
ЭЕ, V ЭЕ
%n ^ ^
I, Auf =---1 Fn +TnWn — явная схема, I — единичная матрица.
Vi ЭЕ,
В неявной схеме можно выделить слагаемые, отличающие ее от явной:
V£ F" + г"№" — £ R nbr AuST-'. (8)
1 dEi Щ________^_________
сглаживатель
Эти дополнительные слагаемые исчезают при сходимости к стационарному решению. Но в процессе сходимости они обеспечивают устойчивость схемы, сглаживая возмущения, порождаемые неустойчивым при больших шагах по времени явным оператором. Далее эти дополнительные слагаемые будут называться неявным сглаживателем.
В общем случае трехмерных уравнений Рейнольдса получающаяся система линейных алгебраических уравнений будет иметь ленточную матрицу с количеством ненулевых диагоналей, соответствующим неявному шаблону схемы. Каждый блок строк матриц соответствует одной ячейке. Граничные условия влияют на строки, соответствующие приграничной ячейке. В общем случае, матрица не удовлетворяет известным достаточным условиям сходимости метода Гаусса — Зейделя [18]. Однако, если применить небольшое количество итераций, как это сделано в [28, 29], можно получить сходимость к стационарному решению. При этом происходит поочередный обход всех ячеек расчетной области и подставление в формулу для текущей ячейки итерационного значения из соседних ячеек. Для учета приращений из соседних ячеек другого блока на каждой итерации выполняется согласование на границах блоков.
Вышеизложенное позволяет хранить все соответствующие ненулевые блоки матрицы привязанными к каждой конкретной ячейке и не хранить общую матрицу.
КОМБИНИРОВАННЫЙ МЕТОД С ЯВНОЙ И НЕЯВНОЙ ЧАСТЬЮ
При форме записи (8) явная схема Годунова — Колгана [24, 38], имеющая второй порядок аппроксимации по пространству, первый порядок по времени и устойчивая при локальном числе Куранта менее 0.5, входит в (8) как упрощенный вариант неявной схемы (схемы со сглаживателем), из которого исключены повышающие устойчивость «неявные» добавки. При локальном выполнении условия устойчивости оказывается возможным не вычислять эти добавки, т. е. матрицу R — I и матрицы R NBR, фактически применяя для таких ячеек явную схему без общего нарушения алгоритма.
Обычно структурированная расчетная сетка имеет сгущения к областям, в которых предполагается протекание наиболее интересных физических явлений. Например, у расчетной сетки для аэродинамического профиля это положение точки растекания на носике, разрешение скачка уплотнения и разрешение схода вихревой пелены на конце профиля.
При изучении вязких течений моделирование пограничного слоя требует многократного сгущения расчетной сетки к поверхностям с установленными на них условиями прилипания. Такие «погранслойные» ячейки отличаются размером от крупных ячеек вдали от тела, на много порядков. Для встраивания явной схемы Годунова — Колгана — Родионова [24, 38, 39], имеющей второй порядок аппроксимации и по пространству, и по времени и устойчивой при локальном числе Куранта менее 1, требуется модификация суммы потоков, входящей в правую часть СЛАУ — локальное введение процедуры «предиктор».
Таким образом, неявный сглаживатель, повышающий устойчивость, но повышающий и диссипативность схемы, и требующий большего количества вычислений на совершение одного шага по времени, применяется лишь там, где схема без него неустойчива. Такой подход вбирает в себя лучшие качества явной и неявной схем.
Для ускорения сходимости в областях, где применяется явная схема, можно использовать технологию локального шага по времени. Это позволит возмущениям быстрее доходить до границ расчетной области по крупным ячейкам расчетной сетки вдали от обтекаемой поверхности.
Кроме того, для еще большего ускорения в буферных областях, интересующих нас только с точки зрения невозвращения возмущений, применяется неявная схема с локальным шагом по времени.
( \ I + Л—I
V сглаживатель Jj
I + R — I Auf =-
Введем определения: назовем «пристеночными» такие ячейки, в которых величина заданного глобального шага по времени превышает локальное условие устойчивости для явной схемы.
Назовем «основными» такие ячейки, в которых величина заданного глобального шага по времени не превышает локальное условие устойчивости для явной схемы и при этом превышает 0.01 (одну сотую) от последнего.
Оставшиеся ячейки (в которых величина заданного глобального шага по времени не превышает 0.01 величины, обеспечивающей локальное условие устойчивости для явной схемы) назовем «буферными».
Единственным параметром, который управляет величиной шага по времени, является заданная величина глобального шага по времени. Алгоритмы выбора величины этого параметра для каждого конкретного типа задач вычислительной аэродинамики выходят за рамки данной статьи.
ТЕСТИРОВАНИЕ МЕТОДА
Описанный выше метод был реализован в виде программы-солвера COMGLEI (Combination of Global and Local tau type with Explicit and Implicit schemes). Эта программа является модификацией неявного блока программы ZEUS, разработанного совместно автором настоящей статьи и С. В. Михайловым. После модификаций удалось добиться совпадения результатов программ COMGLEI и V3Solver.
Тест 1 — турбулентный пограничный слой (ПС) на пластине. Для теста выбран режим M = 0.8. Число Рейнольдса, соответствующее длине пластины, ReL = 22.8 106. Поверхность пластины считалась теплоизолированной. Была построена расчетная сетка, содержащая 165 ячеек вдоль пластины, 200 ячеек поперек пластины, из них 96 ячеек — внутри пограничного слоя
(в конце пластины), при этом размер пристеночных ячеек обеспечивает y + е^10-4;10-3 ), где
y + = y/l*, u + = u/v*, v* = yj\rw\/pw , l* =Vw/v* — универсальные параметры ПС [42] вычисляются на основе данных из расчета. Для оценки точности полученных результатов используется теоретическая зависимость, описанная в [42, 43]:
w+= y+, 0 < y + < 7
y+ = u ++ e-2 \
exp
(.4u +)-1 - 0.4w + -
(.4u +)2 (.4u +)3 (.4u +)4
u =----ln(y ) + 5 +
0.4
П ( x) 0.4
2sin
w
2 V* ô
y + > 30.
16.16
k 7 < y+< 30,
Сравнение результатов (рис. 1) показывает, что максимальное отличие местной величины скорости решений в профиле пограничного слоя в конце пластины, полученных расчетом по явной схеме с локальным шагом по времени (солвером V3Solver) и по описанной выше технологии (солвером СОМОЬБ1), составляет менее 0.004 относительно величины местной скорости и находится в глубине пограничного слоя.
Это различие результатов можно объяснить тем, что при получении решения по явной схеме при стремлении к сходимости медленно сокращается длина ламинарного участка. И, несмотря на то, что расчет по явной схеме продолжался в 5 раз дольше, сходимость решения все еще не достигнута. При использовании СОМОЬБ1 сходимость по интегральному коэффициенту силы трения с^ с точностью 1% достигается в 27 раз быстрее (рис. 2), чем при использовании явной
схемы с локальным шагом по времени. На оси абсцисс отложено время проведения расчета, измеряемое в условных единицах, равных среднему времени расчета одного шага по времени солвером V3Solver.
0.1 1 ю 100 1000 -v+ Условное процессорное время
Рис. 1. Профиль безразмерной скорости в сечении Рис. 2. Сравнение сходимости коэффициента силы трех = 1200 мм (тест 1) ния (тест 1)
Тест 2 — профиль NACA 0012. Для сравнения с расчетом используются экспериментальные данные [44].
Длина хорды профиля L = 1 м. В эксперименте имелся ламинарно-турбулентный переход на расстоянии 0.05L. Расчетная сетка (рис. 3) состоит из четырех блоков и имеет сгущения к поверхности профиля, размер пристеночных ячеек обеспечивает у+ е(0.02; 0.7). Поверхностная
сетка имеет сгущения к носку и к задней кромке профиля. Поверхность моделируется граничным условием «теплоизолированная поверхность с прилипанием». На внешних границах расчетной области — мягкое граничное условие набегающего потока, основанное на инвариантах Римана.
Для теста выбран режим: M = 0.799, а = 2.26°, Re = 9 -106. Турбулентность в точке, удаленной от границ расчетной области на такое же расстояние, что и профиль, q = 0.92 м/с, w = 0.77 Гц при скорости набегающего потока примерно 270 м/с. На этом режиме обтекание носит отрывный характер. Следует отметить, что модель турбулентности q-w не позволяет качественно моделировать отрыв.
Примерно на середине верхней поверхности профиля расположен скачок уплотнения. Из-под него начинается область отрывного течения. На рис. 4 приведено сравнение распределения коэффициента давления, полученного в расчетах как по явной схеме с локальным шагом по времени, так и по рассматриваемой технологии. Кроме того, приведены экспериментальные данные. Различия между расчетами достигают величины 0.001 по коэффициенту давления ср. На рис. 5
приведена история сходимости коэффициента подъемной силы в расчетах. При использовании COMGLEI на рассматриваемом отрывном режиме сходимость по коэффициенту подъемной силы Су с точностью 0.01 достигается в 3 раза быстрее, чем при использовании явной схемы
Рис. 3. Расчетная сетка около профиля (тест 2)
Рис. 4. Распределение коэффициента давления на профиле при М = 0.799, а = 2.26° (тест 2)
с локальным шагом по времени. Однако во всех расчетах имеются незатухающие колебания с амплитудой менее 0.01.
Также были проведены расчеты для безотрывного режима М = 0.8, а = 0, Яе = 9 106. Для этого режима можно сравнить с экспериментом только величину коэффициента сопротивления. Экспериментальные данные разбросаны в диапазоне [0.01, 0.017]. И солвер У38о1уег, и солвер СОМОЬБ1 дают величину 0.0166. На рис. 6 приведены графики сходимости. При использовании СОМОЬБ1 на рассматриваемом безотрывном режиме сходимость по су с точностью 0.001 дости-
Рис. 5. Сравнение сходимости коэффициента подъемной силы профиля при М = 0.799, а = 2.26° (тест 2)
гается в 20 раз быстрее, чем при использовании явной схемы с локальным шагом по времени.
Предложенная в данной работе комбиниро- Рис. 6. Сравнение сходимости коэффициента подъем-ванная технология в тестах 1 и 2 совпадает с неяв- ной силы профиля при М = 0.8, а = 0 (тест 2) ной схемой с глобальным шагом по времени, так
как заданная величина глобального шага по времени была столь большой, что во всей расчетной области не было ячеек, для которых выполнялось бы локальное условие устойчивости явной схемы. Поэтому сравнения с отдельной реализацией неявной схемой с глобальным шагом по времени не проводилось.
Тест 3 — компоновка «фюзеляж — крыло». Многоблочная структурированная расчетная сетка взята на сайте http://aaac.larc.nasa.gov/tsab/cfdlarc/aiaa-dpw/Workshop4/workshop4.html и содержит примерно 3 500 000 ячеек. Сетка была модифицирована — изменено сгущение к поверхности с прилипанием, а также перестроена поверхностная сетка на крыле и горизонтальном оперении. Удаленная часть сетки прорежена и отделена условием несогласованной стыковки блоков «Connect» [34].
Сетка отмасштабирована к размерам экспериментальной модели. Математическая модель представляет собой правую половину ЛА. Геометрические параметры для получения аэродинамических коэффициентов и сечения сравнения локальных параметров взяты из статьи [45]. Поверхностная расчетная сетка представлена на рис. 7.
Рис. 7. Поверхностная расчетная сетка (тест 3)
Для сравнения расчетных и экспериментальных данных взят экспериментальный пуск 106 (число M = 0.85, полное давление Ро = 301 800 Па, полная температура То = 142.5 K,
Re = 19.8 106, скоростной напор q^ = 95 185 Па). Турбулентность на границе расчетной области
(q = 1.7 м/с, W = 50 Гц) была задана таким образом, что в районе модели величина q составляет около 0.3% от скорости набегающего потока.
На рис. 8 приведено сравнение интегральных характеристик, полученных в расчетах как по явной схеме с локальным шагом по времени, так и по рассматриваемой технологии. Кроме того, приведены экспериментальные данные. Различия между расчетами достигают 0.0001 по коэффициенту сопротивления cx. Отличие коэффициента cx, полученного в расчете, от экспериментального значения составляет 0.01.
На рис. 9 представлены графики сходимости расчетов по трем методам: по явной схеме с локальным шагом по времени (Local), по неявной схеме с глобальным шагом по времени (global implicit) и по предложенному комбинированному методу (COMGLEI). Заданная величина глобального шага по времени превышала условие устойчивости для явной схемы в 105 раз и в расчете «global implicit», и в расчете «COMGLEI».
При использовании COMGLEI на рассматриваемом режиме сходимость по су с точностью
0.01 достигается в 5 раз быстрее, чем при использовании явной схемы с локальным шагом по времени, с точностью 0.001 — в 2 раза быстрее. Скорость сходимости в расчетах по неявной схеме с глобальным шагом по времени уступает скорости сходимости в расчетах по явной схеме с локальным шагом по времени.
Рис. 8. Сравнение интегральных характеристик, Рис. 9. Сравнение сходимости по коэффициенту подъ-полученных расчетами и экспериментом (тест 3) емной силы (тест 3)
Рис. 10. Зоны применения различных схем в расчетах теста 3 методом
COMGLEI
На рис. 10 изображено распределение зон, в которых работают различные типы схем, применяемых в методе COMGLEI. Для простоты визуализации приведена плоскость симметрии. В середине находится модель исследуемого ЛА. Слой вокруг модели — зона применения неявной численной схемы с глобальным шагом по времени. Далее — зона применения явной схемы с локальным шагом по времени. На периферии и в вытянутой области, расположенной на подветренной стороне (окруженной зоной применения явной схемы), — зона применения неявной схемы с локальным шагом по времени. Появление упомянутой вытянутой области связано с особенностями структурированной расчетной сетки.
ЗАКЛЮЧЕНИЕ
Предложен комбинированный метод, основанный на схеме Годунова — Колгана — Родионова. Метод сочетает локальное использование неявной схемы с заданным глобальным шагом по времени в окрестности поверхностей с условием прилипания, явной схемы с локальным шагом по времени в областях, где заданный глобальный шаг не превосходит условие устойчивости явной схемы, и неявной схемы с локальным шагом по времени в буферных областях, удаленных от изучаемого тела.
Предложенный комбинированный метод позволил ускорить получение стационарных решений уравнений Рейнольдса по сравнению c имевшимися в ППП EWT-ЦАГИ явными методами, а именно:
в тестовой задаче «пограничный слой на пластине» — в 27 раз (сходимость по с^- с точностью 1%);
в тестовой задаче «профиль NACA0012» на отрывном режиме — в 3 раза (сходимость по су с точностью 0.01), на безотрывном режиме — в 20 раз (сходимость по су с точностью 0.001);
в тестовой задаче «компоновка фюзеляж — крыло» — в 5 раз (сходимость по су с точностью 0.01), либо в 2 раза (сходимость по су с точностью 0.001);
Перечисленные результаты по ускорению расчета были получены на кластере «LEGO» НИО-1 ЦАГИ. Характеристики кластера детально описаны в [46].
1.Neyland V. Ya., Bosniakov S. M., Glazkov S. A., Ivanov A., Mat-yash S. V.Mikhailov S. V.Vlasenko V. V. Conception of electronic wind tunnel and first results of its implementation // Progress in Aerspace Sciences. 2001. V. 37, p. 121 —145.
2. Босняков С. М., Власенко В. В., Курсаков И. А., Михайлов С. В., Квест Ю. Задача интерференции оживального тела вращения с державкой аэродинамической трубы и особенности ее решения с использованием ЭВМ // Ученые записки ЦАГИ. 2011. Т. XLII, № 3, с. 25—40.
3. Сб. статей «Практические аспекты решения задач внешней аэродинамики двигателей летательных аппаратов в рамках осредненных по времени уравнений Навье — Стокса» // Труды ЦАГИ. 2007, вып. 2671.
4. Енгулатова М. Ф., Матяш С. В. Пакет EWT — программная реализация технологии аэродинамического эксперимента // Материалы XVIII школы-семинара ЦАГИ «Аэродинамика летательных аппаратов». — 2007.
5. Власенко В. В., Михайлов С. В. Программа ZEUS для расчета нестационарных течений в рамках подходов RANS и // Материалы XVIII школы-семинара ЦАГИ «Аэродинамика летательных аппаратов». — 2009.
6. Blazek J. Computational fluid dynamics: principles and applications. — New York: Elsevier, 2001.
7. Федоренко Р. П. Многосеточный метод для схем конечного элемента // ЖВМ и МФ. 1961. Т. 1, № 5.
8. Trottenberg U., Oosterlee C., Schuller A. Multigrid. — Academic Press, MPG Books, Cornwall, 2001.
9.Yoon S. & Kwa k D. Implicit methods for the Navier — Stokes equations // Computing Systems in Engineering, 1990. V. 1, N 2—4, p. 535—547.
10. Venkatakrishnan V. Newton Solution of inviscid and viscous problems // AIAA J. 1989. V. 27, p. 885—891.
11. Hestenes M. R., Stiefel E. L. Methods of conjugate gradients for solving linear systems // J. Res. Nat. Bur. Stand. 1952. V. 49, p. 409.
12. Sonneveld P. CGS, a fast Lanczos-type solver for nonsymmetric linear systems // SIAM J. Scientific Statistsics and Computing. 1989. V. 10, p. 36—52.
13. Sad Y., Schulz M. H. GMRES: a generalized minimum residual algorithm for solving nonsymmetric linear systems // SIAM J. Scientific and Statistical Computing. 1986. V. 7, p. 856—869.
14. Zingg D.Pueyo A. An eficient Newton-GMRES solver for aerodynamic computations // AIAA Paper 97-1955. 1997.
15. Venkatakrishnan V. Preconditioned conjugate gradient methods for the compressible Navier — Stokes equations // AIAA J. 1991. V. 29, p. 1092—1110.
16. Beam R. M., Warming R. F. An implicit finite-difference algorithm for hyperbolic systems in conservation-law form // J. of Computational Physics. 1976. V. 22, p.87—110.
17. S t e g e r J. L., Warming R. F. Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods // J. of Computational Physics. 1981. V. 40, p. 263—293.
18. Ortega M. J. Introduction to parallel and vector solution of linear systems. — New York: Plenum Press, 1988.
19. Yoon S., Jameson A. An LU-SSOR scheme for the Euler and Navier — Stokes equations // AIAA J. 1988. V. 26, p. 1025—1026.
20. B a t i n a J. T. Implicit upwind solution algorithms for three-dimensional unstructured meshes // AIAA J. 1993. V. 31, p. 801—805.
21. Anderson W. K, Bonhaus D. L. An implicit upwind algorithm for computing turbulent flows on unstructured grids // Computers & Fluids. 1994. V. 23, p. 1 —21.
22. Ferziger J., Peric M. Computational methods for fluid dynamics. 3rd rev.ed. — New York: Springer, 2002.
23. Годунов С. К. Разностный метод численного расчета разрывных решений гидродинамики // Математический сборник. 1959. Вып. 47(89), № 3, с. 271—306.
24. Годунов С. К., Забродин А. В., Иванов М. Я., Крайко А. Н., Прокопов Г. П. Численное решение многомерных задач газовой динамики. — М.: Наука, 1976.
25. Иванов М. Я., Нигматулин Р. З. Неявная схема С. К. Годунова повышенной точности для численного интегрирования уравнений Эйлера // ЖВМ и МФ. 1987. Т. 27, № 11, с. 1725 — 1735.
26. Иванов М. Я., Крупа В. Г., Нигматулин Р. З. Неявная схема С.К. Годунова повышенной точности для интегрирования уравнений Навье — Стокса// ЖВМ и МФ. 1989. Т. 29, № 6, с. 888—901.
27. R o e P. L. Approximate Riemann solvers, parameter vectors, and difference schemes // J. of Comput. Phys. 1981. V. 43, p. 357—372.
28. Топеха Е. А., Копчёнов В. И. Неявная релаксационная конечно-разностная схема для системы уравнений Эйлера // Методы исследования гиперзвуковых летательных аппаратов. Ежегодная научная школа-семинар ЦАГИ «Механика жидкости и газа». — Сб. докладов. — Изд. ЦАГИ. 25 февраля — 1 марта 1992. Часть 3.
29. Топеха Е. А., Копчёнов В. И. Неявная релаксационная конечно-разностная схема для системы уравнений Навье — Стокса // Методы исследования гиперзвуковых летательных аппаратов. Ежегодная научная школа-семинар ЦАГИ «Механика жидкости и газа». — Сб. докладов. — Изд. ЦАГИ. 1994, с. 9.1—9.10.
30. Программа COBRA. — Свидетельство о регистрации № 2011615671.
31. Wilcox D. C. Progress in hypersonic turbulence modeling // AIAA Paper 91-1785,
1991.
32. Власенко В.В. О математическом подходе и принципах построения численных методологий для пакета прикладных программ EWT ЦАГИ. — В сб.: Практические аспекты решения задач внешней аэродинамики двигателей летательных аппаратов в рамках осредненных по времени уравнений Навье — Стокса // Труды ЦАГИ. 2007, вып. 2671, с. 20—85.
33. Spalart P. R.Allmaras S. R. A One-equation turbulence model for aerodynamic flows // AIAA Paper 92-439, 1992.
34. Glowinski R., Golub G. H., Meurant G. A., Periaux J. (Eds.) Domain decomposition methods for partial differential equations. — New York: Hardcover, 1988, p. 431, ISBN: 0-89871-220-3.
35. C o akley T. J. Turbulence modeling methods for the compressible Navier — Stokes equations // AIAA Paper 83-1693, 1983.
36. Co akley T., Hsieh T. Comparison between implicit and hybrid methods for the calculation of steady and unsteady inlet flows // AIAA Paper 85-1125, 1985.
37. Михайлов С. В. Адаптация коэффициентов q-ю модели турбулентности к особенностям течения // Материалы XX Школы-семинара ЦАГИ «Аэродинамика летательных аппаратов», 2009.
38. Колган В. П. Применение принципа минимальных значений производной к построению конечно-разностных схем для расчета разрывных решений газовой динамики // Ученые записки ЦАГИ. 1972. Т. 3, № 6, с. 68—77.
39. Родионов А. В. Монотонная схема второго порядка аппроксимации для маршевых расчетов неравновесных потоков // ЖВМ и МФ. 1987. Т. 27, № 4.
40. E l i a s s o n P. EDGE: A Navier — Stokes solver for unstructured grids // Proc. Finite Volumes for Complex Applications III. 2002, p. 527—534.
41. Оран Э., Борис Дж. Численное моделирование реагирующих потоков. — М: Мир, 1990.
42. Хлопков Ю. И., Жаров В. К., Горелов А. Л. Когерентные структуры в ТПС. — М: МФТИ, 2002.
43. Vlasenko V., Shiryaeva A. Numerical simulation of non-stationary propagation of combustion along a duct with supersonic flow of a viscid gas // Proceedings of the Institution of Mechanical Engineers, Part G: J. of Aerospace Engineering, 2012.
44. Terry L. Holst viscous transonic airfoil workshop compendium of results //AIAA Paper 87-1460.
45. Rivers M. B., Dittberner A. Experimental Investigation of the NASA Common Research Model in the NASA Langley National Transonic Facility and NASA Ames 11-Ft Transonic Tunnel //AAIA Paper 2011-1126.
46. Кажан Е. В., Курсаков И. А., Лысенков А. Л. Установка, наладка и пуск в эксплуатацию новой вычислительной системы — кластер «LEGO» // Материалы XX школы-семинара ЦАГИ «Аэродинамика летательных аппаратов», 2009.
Рукопись поступила 6/II2012 г.
